Chapter 9 — Deconvolution & CLEAN¶

!!! info "Before you start" Prerequisites: Ch 8 (Aperture Synthesis & the uv-plane) · Maths Lab: Lab A (Fourier & Convolution) · ~50 min · Intermediate

In Chapter 8 we learned that an interferometer is a Fourier machine: each baseline measures one visibility $V(u, v)$ — one Fourier component of the sky — and a real array samples that Fourier plane only where it has baselines. The image we recovered, the dirty image, was the true sky $I(l,m)$ convolved with the dirty beam $B(l,m)$, the point-spread function set by the holes in our $uv$ coverage:

$$ I_\mathrm{D}(l, m) \;=\; I(l, m)\;\ast\;B(l, m). $$

That dirty image was riddled with sidelobes — rings and spokes thrown across the field by every bright source. It is not yet a science-quality picture. This chapter is about undoing that convolution: deconvolution, and specifically the algorithm that made aperture-synthesis imaging practical, Högbom's CLEAN (1974).

Learning goals¶

By the end of this chapter you will be able to:

• State the deconvolution problem: given the dirty image $I_\mathrm{D}$ and the dirty beam $B$, recover an estimate of the true sky $I$.
• Explain why naive Fourier division fails (division by zero in the unsampled holes) and why we need a constrained, iterative method instead.
• Walk through Högbom's CLEAN step by step: find the peak, subtract a scaled dirty beam, accumulate a clean component, repeat — then restore the components with a clean (well-behaved) beam.
• Run jansky.interferometry.hogbom_clean on a real dirty image and watch the sources emerge from the sidelobes.
• See, qualitatively, how CASA's tclean generalises this idea (a forward pointer to Chapter 12).

The history: Högbom and CLEAN¶

By the early 1970s, Earth-rotation synthesis arrays could fill the $uv$-plane well enough to make real maps — but the maps were dirty. A bright source spread a corona of sidelobes across the whole field, burying fainter sources and making the images hard to trust. The sampling function $S(u,v)$ was full of holes, so the synthesised beam had nasty sidelobes, and there was no obvious way to divide them back out.

Jan Högbom, working at Cambridge and Westerbork, proposed a beautifully simple, almost greedy, idea. His paper is:

Högbom, J. A. (1974). Aperture Synthesis with a Non-Regular Distribution of Interferometer Baselines. Astronomy & Astrophysics Supplement Series 15, 417. ADS

The insight: model the sky as a collection of point sources. If the true sky were a single point, the dirty image would be exactly one copy of the dirty beam, centred on that point. So: look at the brightest pixel in the dirty image, assume some of that flux is a real point source, record it, and subtract the dirty beam (scaled and shifted to that pixel) from the image. That removes not just the peak but all of its sidelobes too. Repeat on what's left. Each iteration peels off one more clean component, and the leftover residual gets flatter and flatter until it is just noise.

Because most radio sources really are dominated by compact structure, this "deconvolve-by-modelling-as-points" trick works astonishingly well. CLEAN turned synthesis imaging from a curiosity into the workhorse of radio astronomy. Every modern imager — including CASA's tclean, which we will meet in Chapter 12 — is a descendant of these few lines from 1974.

The standard modern treatment is Thompson, Moran & Swenson, Interferometry and Synthesis in Radio Astronomy (3rd ed., open access), and the free Condon & Ransom, Essential Radio Astronomy (online), both of which devote a chapter to deconvolution.

The deconvolution problem¶

We are handed two arrays and want a third. We have the dirty image $I_\mathrm{D}$ and the dirty beam $B$, and we know they are related to the unknown true sky $I$ by a convolution:

$$ I_\mathrm{D} \;=\; I \;\ast\; B . $$

In Fourier space convolution becomes multiplication, so naively you might try to recover $I$ by dividing:

$$ \hat I \;=\; \mathcal{F}^{-1}\!\left[\frac{\mathcal{F}[I_\mathrm{D}]} {\mathcal{F}[B]}\right] . $$

This fails. The transform of the beam, $\mathcal{F}[B] = S(u,v)$, is the sampling function — and it is zero everywhere the array had no baseline. You would be dividing by zero across every hole in the $uv$-plane. Even a tiny bit of noise gets amplified without bound. The unmeasured visibilities are simply missing information, and no amount of arithmetic conjures them back.

So deconvolution here is ill-posed: many different skies are consistent with the data. We need a constraint — some prior belief about what the sky looks like — to pick a sensible one. CLEAN's constraint is: the sky is a modest number of point sources.

Högbom's CLEAN, step by step¶

Given the dirty image and dirty beam, CLEAN builds up a model (a list of clean components) like this:

1. Find the peak. Locate the pixel $(i, j)$ of largest absolute value in the current residual image (which starts equal to the dirty image).
2. Record a component. Add a fraction $g$ of that peak — the loop gain, typically $g \approx 0.1$ — as a point source in the model at $(i, j)$: $\;\Delta = g \cdot I_\mathrm{res}(i, j)$.
3. Subtract its beam. Subtract a copy of the dirty beam, scaled by $\Delta$ and shifted so its peak sits on $(i, j)$, from the residual. This removes the component and all of its sidelobes.
4. Repeat from step 1 until you hit a maximum number of iterations or the residual peak drops below a threshold (ideally, the noise level).

The small loop gain is what keeps CLEAN stable: instead of claiming all of a peak at once (which over-fits and creates artefacts), it nibbles away a little at a time, letting overlapping sidelobes from many components settle out gradually.

This is exactly what jansky.interferometry.hogbom_clean implements — open src/jansky/interferometry.py and you will see these four steps in about a dozen lines of NumPy.

Restoring: the clean beam¶

After CLEAN converges you hold two things: a model of delta-function components, and a residual image of leftover noise and un-modelled flux. You could just show the model — but a set of mathematical delta functions is visually misleading: it claims infinite resolution your data never had, and any small position errors in the components look jarring.

So we restore. We convolve the clean-component model with a clean beam — a smooth Gaussian fitted to the central lobe of the dirty beam, with the sidelobes thrown away — and then add the residual back on top:

$$ I_\mathrm{clean} \;=\; \big(\,\text{model}\,\ast\,B_\mathrm{clean}\,\big) \;+\; I_\mathrm{res} . $$

The result is an honest image: it has the resolution the array actually achieved (set by the main lobe), the sidelobes are gone, and adding the residual back preserves any flux CLEAN did not get around to modelling (and lets you see left-over noise). This restored image is what gets published.

Setup¶

We reuse the same helpers as Chapter 8 — jansky.interferometry for the synthesis machinery and jansky.plotting for consistent figures. The only new star of the show is interferometry.hogbom_clean.

In [ ]:

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter

from jansky import interferometry, plotting
from jansky.interferometry import (
    uv_coverage,
    grid_visibilities,
    dirty_beam,
    dirty_image,
    hogbom_clean,
    CleanResult,
)
from jansky.plotting import use_jansky_style, plot_uv_coverage, show_image

use_jansky_style()
rng = np.random.default_rng(9)  # reproducible 'randomness' for this chapter

print("CLEAN helper:", hogbom_clean.__name__)
print("Returns a    :", CleanResult.__name__, "(model, residual, components)")

import numpy as np import matplotlib.pyplot as plt from scipy.ndimage import gaussian_filter from jansky import interferometry, plotting from jansky.interferometry import ( uv_coverage, grid_visibilities, dirty_beam, dirty_image, hogbom_clean, CleanResult, ) from jansky.plotting import use_jansky_style, plot_uv_coverage, show_image use_jansky_style() rng = np.random.default_rng(9) # reproducible 'randomness' for this chapter print("CLEAN helper:", hogbom_clean.__name__) print("Returns a :", CleanResult.__name__, "(model, residual, components)")

A model sky and its dirty image (recap of Chapter 8)¶

We rebuild the same ingredients we used in Chapter 8, so this chapter stands on its own: a small, deliberately irregular array, swept through Earth rotation to fill the $uv$-plane; a model sky of two point sources plus one soft extended blob; and the dirty beam and dirty image that result.

In [ ]:

# A small, deliberately irregular 8-antenna array (positions in metres).
antenna_xy = np.array(
    [
        [0.0, 0.0],
        [140.0, 20.0],
        [-90.0, 130.0],
        [210.0, -110.0],
        [-180.0, -60.0],
        [60.0, 240.0],
        [-250.0, 150.0],
        [300.0, 90.0],
    ]
)

# Earth-rotation synthesis: track the source across +/- 6 h of hour angle,
# near the north celestial pole, to fill the uv-plane.
declination = np.deg2rad(80.0)
hour_angles = np.deg2rad(15.0 * np.linspace(-6.0, 6.0, 120))
uv_track = uv_coverage(antenna_xy, hour_angles=hour_angles, declination=declination)

npix = 256


def gaussian_blob(npix, row, col, amp, sigma):
    # A round Gaussian source centred at (row, col) on a npix x npix grid.
    yy, xx = np.mgrid[0:npix, 0:npix]
    r2 = (yy - row) ** 2 + (xx - col) ** 2
    return amp * np.exp(-r2 / (2.0 * sigma**2))


# True sky: two point sources + one extended blob.
sky = np.zeros((npix, npix))
sky[120, 110] = 1.0  # bright point source
sky[150, 165] = 0.6  # fainter point source
sky += gaussian_blob(npix, row=110, col=150, amp=0.4, sigma=6.0)  # extended blob

# Grid the uv samples, then form the dirty beam (PSF) and dirty image.
sampling = grid_visibilities(uv_track, npix=npix)
beam = dirty_beam(sampling)
dirty = dirty_image(sky, sampling)

print("sky    peak =", round(float(sky.max()), 3))
print("beam   peak =", round(float(beam.max()), 3), "(normalised)")
print("dirty  peak =", round(float(dirty.max()), 3))

# A small, deliberately irregular 8-antenna array (positions in metres). antenna_xy = np.array( [ [0.0, 0.0], [140.0, 20.0], [-90.0, 130.0], [210.0, -110.0], [-180.0, -60.0], [60.0, 240.0], [-250.0, 150.0], [300.0, 90.0], ] ) # Earth-rotation synthesis: track the source across +/- 6 h of hour angle, # near the north celestial pole, to fill the uv-plane. declination = np.deg2rad(80.0) hour_angles = np.deg2rad(15.0 * np.linspace(-6.0, 6.0, 120)) uv_track = uv_coverage(antenna_xy, hour_angles=hour_angles, declination=declination) npix = 256 def gaussian_blob(npix, row, col, amp, sigma): # A round Gaussian source centred at (row, col) on a npix x npix grid. yy, xx = np.mgrid[0:npix, 0:npix] r2 = (yy - row) ** 2 + (xx - col) ** 2 return amp * np.exp(-r2 / (2.0 * sigma**2)) # True sky: two point sources + one extended blob. sky = np.zeros((npix, npix)) sky[120, 110] = 1.0 # bright point source sky[150, 165] = 0.6 # fainter point source sky += gaussian_blob(npix, row=110, col=150, amp=0.4, sigma=6.0) # extended blob # Grid the uv samples, then form the dirty beam (PSF) and dirty image. sampling = grid_visibilities(uv_track, npix=npix) beam = dirty_beam(sampling) dirty = dirty_image(sky, sampling) print("sky peak =", round(float(sky.max()), 3)) print("beam peak =", round(float(beam.max()), 3), "(normalised)") print("dirty peak =", round(float(dirty.max()), 3))

In [ ]:

# The starting point: the true sky (which we never observe directly),
# the dirty beam (PSF), and the dirty image CLEAN must deconvolve.
fig, axes = plt.subplots(1, 3, figsize=(15, 5.2))
show_image(sky, ax=axes[0], title="True sky (unknown to the telescope)")
show_image(beam, ax=axes[1], title="Dirty beam  B  (PSF, with sidelobes)")
show_image(dirty, ax=axes[2], title="Dirty image  I_D = sky * B")
plt.tight_layout()
plt.show()

# The starting point: the true sky (which we never observe directly), # the dirty beam (PSF), and the dirty image CLEAN must deconvolve. fig, axes = plt.subplots(1, 3, figsize=(15, 5.2)) show_image(sky, ax=axes[0], title="True sky (unknown to the telescope)") show_image(beam, ax=axes[1], title="Dirty beam B (PSF, with sidelobes)") show_image(dirty, ax=axes[2], title="Dirty image I_D = sky * B") plt.tight_layout() plt.show()

Look at the dirty image: the two point sources are there, but smeared and wrapped in sidelobe rings, and the soft blob is half-buried in the ripples thrown off by the bright source. This is the input to CLEAN.

Running CLEAN¶

We hand the dirty image and dirty beam to hogbom_clean with a conservative loop gain of $g = 0.1$ and a generous iteration budget. It returns a CleanResult with three fields:

• model — the clean components painted onto the pixel grid (delta functions);
• residual — what is left of the dirty image after subtracting every component;
• components — the ordered list of (row, col, flux) that CLEAN found.

In [ ]:

gain = 0.1  # loop gain: subtract 10% of each peak per iteration
n_iter = 500  # maximum number of CLEAN iterations

result = hogbom_clean(dirty, beam, gain=gain, n_iter=n_iter)

print(f"CLEAN found {len(result.components)} components in <= {n_iter} iterations")
print(
    f"residual peak = {np.abs(result.residual).max():.4f}  "
    f"(dirty peak was {np.abs(dirty).max():.4f})"
)
print("brightest few components (row, col, flux):")
for row, col, flux in result.components[:5]:
    print(f"  ({row:3d}, {col:3d})  flux = {flux:.4f}")

gain = 0.1 # loop gain: subtract 10% of each peak per iteration n_iter = 500 # maximum number of CLEAN iterations result = hogbom_clean(dirty, beam, gain=gain, n_iter=n_iter) print(f"CLEAN found {len(result.components)} components in <= {n_iter} iterations") print( f"residual peak = {np.abs(result.residual).max():.4f} " f"(dirty peak was {np.abs(dirty).max():.4f})" ) print("brightest few components (row, col, flux):") for row, col, flux in result.components[:5]: print(f" ({row:3d}, {col:3d}) flux = {flux:.4f}")

The very first components should land right on top of our two input point sources at (120, 110) and (150, 165) — CLEAN has found the sources. Notice that no single component carries the full source flux: because the loop gain is 0.1, the bright source is rebuilt from many small components stacked at (and near) the same pixel. That is CLEAN working as intended.

Restoring with a clean beam¶

The raw model is a field of delta functions — visually harsh and overclaiming resolution. We restore it by convolving with a clean beam: a smooth Gaussian matched to the central lobe of the dirty beam.

A full imager (CASA) fits an elliptical Gaussian to the dirty beam's main lobe. For this teaching demo a round Gaussian is plenty — we estimate its width from the half-width of the dirty beam's central peak and convolve with scipy.ndimage.gaussian_filter. Then we add the residual back to get the final restored image.

In [ ]:

def clean_beam_sigma(beam):
    # Estimate a round clean-beam sigma from the dirty beam's central lobe.
    # Walk out from the beam centre along a row until the beam first drops to
    # half its peak; that half-width-at-half-maximum (HWHM) fixes the Gaussian
    # sigma via FWHM = 2.355 * sigma.
    c = beam.shape[0] // 2
    row = beam[c, c:]  # half-slice from the peak outward
    half = np.argmax(row < 0.5 * row[0])  # first pixel below half-max
    hwhm = max(half, 1)
    fwhm = 2.0 * hwhm
    return fwhm / 2.3548


sigma_clean = clean_beam_sigma(beam)
print(f"clean-beam sigma ~ {sigma_clean:.2f} px  (FWHM ~ {2.3548 * sigma_clean:.1f} px)")

# Restore: convolve the clean-component model with the clean beam, add residual.
restored_model = gaussian_filter(result.model, sigma=sigma_clean)
restored = restored_model + result.residual

print("restored image peak =", round(float(restored.max()), 3))

def clean_beam_sigma(beam): # Estimate a round clean-beam sigma from the dirty beam's central lobe. # Walk out from the beam centre along a row until the beam first drops to # half its peak; that half-width-at-half-maximum (HWHM) fixes the Gaussian # sigma via FWHM = 2.355 * sigma. c = beam.shape[0] // 2 row = beam[c, c:] # half-slice from the peak outward half = np.argmax(row < 0.5 * row[0]) # first pixel below half-max hwhm = max(half, 1) fwhm = 2.0 * hwhm return fwhm / 2.3548 sigma_clean = clean_beam_sigma(beam) print(f"clean-beam sigma ~ {sigma_clean:.2f} px (FWHM ~ {2.3548 * sigma_clean:.1f} px)") # Restore: convolve the clean-component model with the clean beam, add residual. restored_model = gaussian_filter(result.model, sigma=sigma_clean) restored = restored_model + result.residual print("restored image peak =", round(float(restored.max()), 3))

The CLEAN sequence, in four panels¶

This is the heart of the chapter. From the dirty image, CLEAN extracts a set of clean components, leaving a flat residual, and we restore the components with the clean beam. Compare the restored image with the true sky above: the point sources are sharp and sidelobe-free, and the extended blob is recovered (a little roughly — CLEAN models extended emission as a cluster of points, its one real weakness).

In [ ]:

fig, axes = plt.subplots(2, 2, figsize=(12, 11))

show_image(dirty, ax=axes[0, 0], title="1. Dirty image  (input)")

# Clean components: mark each found component on a blank grid.
comp_rows = [r for r, c, f in result.components]
comp_cols = [c for r, c, f in result.components]
show_image(
    result.model, ax=axes[0, 1], title=f"2. Clean components  ({len(result.components)} found)"
)
axes[0, 1].scatter(
    comp_cols, comp_rows, s=10, facecolors="none", edgecolors="cyan", linewidths=0.6, alpha=0.5
)

show_image(result.residual, ax=axes[1, 0], title="3. Residual  (after subtraction)")
show_image(restored, ax=axes[1, 1], title="4. Restored image  (components * clean beam + residual)")

plt.tight_layout()
plt.show()

fig, axes = plt.subplots(2, 2, figsize=(12, 11)) show_image(dirty, ax=axes[0, 0], title="1. Dirty image (input)") # Clean components: mark each found component on a blank grid. comp_rows = [r for r, c, f in result.components] comp_cols = [c for r, c, f in result.components] show_image( result.model, ax=axes[0, 1], title=f"2. Clean components ({len(result.components)} found)" ) axes[0, 1].scatter( comp_cols, comp_rows, s=10, facecolors="none", edgecolors="cyan", linewidths=0.6, alpha=0.5 ) show_image(result.residual, ax=axes[1, 0], title="3. Residual (after subtraction)") show_image(restored, ax=axes[1, 1], title="4. Restored image (components * clean beam + residual)") plt.tight_layout() plt.show()

The transformation from panel 1 to panel 4 is what CLEAN buys you, and what made synthesis imaging practical. The sidelobes that smeared flux across the dirty image are gone; the residual (panel 3) is flat and structureless, which is the sign of a good CLEAN — it means almost all the real structure has been pulled into the model. A residual still full of sidelobe rings would tell you to keep iterating (or lower the threshold).

Did CLEAN actually recover the sources?¶

Let's make the success quantitative. We sum the clean-component flux landing on (and immediately around) each input point source, and compare the positions of the brightest components with where we put the sources.

In [ ]:

true_sources = {"bright (1.0 Jy)": (120, 110), "faint (0.6 Jy)": (150, 165)}

print("Recovered clean-component flux near each input source:")
for name, (r0, c0) in true_sources.items():
    flux = sum(f for r, c, f in result.components if abs(r - r0) <= 2 and abs(c - c0) <= 2)
    print(f"  {name:16s} at ({r0}, {c0}): summed component flux = {flux:.3f}")

# The single brightest component should sit on the brightest true source.
peak_r, peak_c, _ = max(result.components, key=lambda t: t[2])
print(f"\nBrightest component is at ({peak_r}, {peak_c}); true bright source is at (120, 110).")

true_sources = {"bright (1.0 Jy)": (120, 110), "faint (0.6 Jy)": (150, 165)} print("Recovered clean-component flux near each input source:") for name, (r0, c0) in true_sources.items(): flux = sum(f for r, c, f in result.components if abs(r - r0) <= 2 and abs(c - c0) <= 2) print(f" {name:16s} at ({r0}, {c0}): summed component flux = {flux:.3f}") # The single brightest component should sit on the brightest true source. peak_r, peak_c, _ = max(result.components, key=lambda t: t[2]) print(f"\nBrightest component is at ({peak_r}, {peak_c}); true bright source is at (120, 110).")

Two things to notice. First, the positions are exact: the brightest component sits right on the bright source, and the relative recovered fluxes ($\approx 0.087$ vs $\approx 0.052$) track the true brightness ratio ($1.0 : 0.6$) closely — CLEAN nailed where the sources are and how bright they are relative to each other. Second, the absolute component sums are small because our $uv$ sampling is normalised so the dirty image peaks at only $\approx 0.09$ rather than $1.0$; CLEAN recovers the flux that is present in the dirty image, and the residual (now $\sim 10^{-4}$) shows almost nothing is left behind. Calibrating those numbers to physical janskys is a flux-scale bookkeeping step that real pipelines handle and we set aside for this teaching demo.

Try it yourself¶

These two exercises are scaffolded with # TODO cells. Each runs as-is (it just re-uses the baseline result), but the learning happens when you fill in the # TODO lines and compare.

Exercise 1 — Vary the loop gain¶

The loop gain $g$ controls how aggressively CLEAN claims each peak. A small gain is slow but stable; a large gain converges fast but can over-fit and leave artefacts. Run CLEAN at several gains with a fixed iteration budget and compare the residual peak (lower is better) and the recovered flux.

In [ ]:

# TODO: try a range of loop gains and see how fast/cleanly CLEAN converges.
#       Fill in the gain values, then run hogbom_clean for each and record
#       the residual peak. Smaller gain -> more iterations needed; larger gain
#       -> faster but riskier. What gain gives the flattest residual here?

gains_to_try = [0.05, 0.1, 0.3]  # TODO: add e.g. 0.5, 0.8 and compare

print(f"{'gain':>6}  {'residual peak':>14}  {'n components':>12}")
for g in gains_to_try:
    # TODO: call hogbom_clean(dirty, beam, gain=g, n_iter=500) and report
    #       np.abs(res.residual).max() and len(res.components).
    res = hogbom_clean(dirty, beam, gain=g, n_iter=500)  # <-- baseline; tweak n_iter
    print(f"{g:6.2f}  {np.abs(res.residual).max():14.4f}  {len(res.components):12d}")

# TODO: try a range of loop gains and see how fast/cleanly CLEAN converges. # Fill in the gain values, then run hogbom_clean for each and record # the residual peak. Smaller gain -> more iterations needed; larger gain # -> faster but riskier. What gain gives the flattest residual here? gains_to_try = [0.05, 0.1, 0.3] # TODO: add e.g. 0.5, 0.8 and compare print(f"{'gain':>6} {'residual peak':>14} {'n components':>12}") for g in gains_to_try: # TODO: call hogbom_clean(dirty, beam, gain=g, n_iter=500) and report # np.abs(res.residual).max() and len(res.components). res = hogbom_clean(dirty, beam, gain=g, n_iter=500) # <-- baseline; tweak n_iter print(f"{g:6.2f} {np.abs(res.residual).max():14.4f} {len(res.components):12d}")

Solution

Add a couple of larger gains to gains_to_try and call hogbom_clean for each at a fixed n_iter. With this generous 500-iteration budget every gain drives the residual peak down to the noise floor, so the difference is in how few iterations it took to get there — larger gains claim more flux per peak and converge faster (but on noisier data they would over-fit and leave artefacts).

gains_to_try = [0.05, 0.1, 0.3, 0.5, 0.8]

print(f"{'gain':>6}  {'residual peak':>14}  {'n components':>12}")
for g in gains_to_try:
    res = hogbom_clean(dirty, beam, gain=g, n_iter=500)
    print(f"{g:6.2f}  {np.abs(res.residual).max():14.4f}  {len(res.components):12d}")

Expected: all five gains hit the 500-component budget and reach a residual peak of $\sim 10^{-4}$ or below (e.g. $\approx 0.0001$ at $g=0.05$–$0.1$, dropping toward $\approx 0.0000$ by $g=0.3$–$0.8$). The takeaway: with enough iterations the loop gain mainly trades speed of convergence against stability — small $g$ is slow but safe, large $g$ is fast but risks over-fitting real noise.

Exercise 2 — Vary the number of iterations¶

Watch the residual flatten as CLEAN runs longer. Run CLEAN at increasing n_iter and plot (or print) the residual peak against iteration count. Where does it level off? Once the residual peak reaches the image noise, extra iterations only "clean the noise" — a classic mistake that invents spurious sources.

In [ ]:

# TODO: sweep n_iter and watch the residual peak fall, then plateau.
#       Fill in the iteration counts below, run CLEAN for each, and plot
#       residual_peak vs n_iter. Identify roughly where it stops improving.

iters_to_try = [10, 50, 100]  # TODO: extend, e.g. [10, 50, 100, 300, 800, 2000]

residual_peaks = []
for n in iters_to_try:
    # TODO: res = hogbom_clean(dirty, beam, gain=0.1, n_iter=n)
    #       residual_peaks.append(np.abs(res.residual).max())
    res = hogbom_clean(dirty, beam, gain=0.1, n_iter=n)
    residual_peaks.append(float(np.abs(res.residual).max()))

fig, ax = plt.subplots(figsize=(8, 4.5))
ax.plot(iters_to_try, residual_peaks, "o-", color="#1f77b4")
ax.set_xlabel("CLEAN iterations")
ax.set_ylabel("residual peak")
ax.set_title("Residual peak vs CLEAN iterations  (extend the list and watch it plateau)")
plt.show()

# TODO: sweep n_iter and watch the residual peak fall, then plateau. # Fill in the iteration counts below, run CLEAN for each, and plot # residual_peak vs n_iter. Identify roughly where it stops improving. iters_to_try = [10, 50, 100] # TODO: extend, e.g. [10, 50, 100, 300, 800, 2000] residual_peaks = [] for n in iters_to_try: # TODO: res = hogbom_clean(dirty, beam, gain=0.1, n_iter=n) # residual_peaks.append(np.abs(res.residual).max()) res = hogbom_clean(dirty, beam, gain=0.1, n_iter=n) residual_peaks.append(float(np.abs(res.residual).max())) fig, ax = plt.subplots(figsize=(8, 4.5)) ax.plot(iters_to_try, residual_peaks, "o-", color="#1f77b4") ax.set_xlabel("CLEAN iterations") ax.set_ylabel("residual peak") ax.set_title("Residual peak vs CLEAN iterations (extend the list and watch it plateau)") plt.show()

Solution

Extend iters_to_try to span a wide range and run CLEAN once per value at fixed gain=0.1, recording np.abs(res.residual).max() each time. Plotting the residual peak against iteration count shows it falling steeply at first and then flattening once almost all the real flux has been pulled into the model.

iters_to_try = [10, 50, 100, 300, 800, 2000]

residual_peaks = []
for n in iters_to_try:
    res = hogbom_clean(dirty, beam, gain=0.1, n_iter=n)
    residual_peaks.append(float(np.abs(res.residual).max()))

fig, ax = plt.subplots(figsize=(8, 4.5))
ax.plot(iters_to_try, residual_peaks, "o-", color="#1f77b4")
ax.set_xlabel("CLEAN iterations")
ax.set_ylabel("residual peak")
ax.set_yscale("log")  # the fall spans orders of magnitude
ax.set_title("Residual peak vs CLEAN iterations")
plt.show()

Expected: the residual peak drops fast — roughly $0.041 \to 0.005 \to 0.0004$ at 10, 50, 100 iterations — then plateaus around $\sim 10^{-4}$ by ~300 iterations and barely improves out to 2000. The takeaway: once the residual peak reaches the noise floor, extra iterations only "clean the noise" and start inventing spurious components, so a sensible stopping threshold (here, a few hundred iterations) matters more than a bigger budget.

From Högbom to CASA's tclean (a peek at Chapter 12)¶

The CLEAN we just ran is the original 1974 algorithm. Production imaging today runs in CASA (the Common Astronomy Software Applications package) through its tclean task, which generalises Högbom's idea in several important ways:

• Major/minor cycles. Rather than working entirely in the image plane, tclean alternates between minor cycles (Högbom-like component-finding on a residual image) and major cycles that go back to the measured visibilities to compute an accurate residual, avoiding error build-up from the approximate image-plane beam.
• Better sky models. Clark and Cotton–Schwab CLEAN speed up the minor cycle; multi-scale CLEAN models extended emission with blobs of several sizes instead of points (fixing exactly the weakness we saw on our blob); and multi-frequency synthesis fits spectral structure across a wide band.
• Weighting and masks. Natural vs. uniform vs. Briggs robust weighting trade sensitivity against resolution; user-drawn or automatic masks tell CLEAN where it is allowed to place components.
• Gridding for real skies. W-projection and A-projection handle the wide fields and antenna effects that our flat 2-D demo ignores.

But strip all of that away and the beating heart is still the same four steps: find the peak, subtract a scaled beam, accumulate a component, restore. Högbom (1974) in a dozen lines of NumPy is genuinely the engine inside a tool that images data from the VLA, ALMA and the SKA. We will run tclean on a real measurement set in Chapter 12.

Recap¶

• The dirty image is the true sky convolved with the dirty beam: $I_\mathrm{D} = I \ast B$. Recovering $I$ is the deconvolution problem.
• Naive Fourier division fails because $\mathcal{F}[B] = S(u,v)$ is zero in every unsampled hole — the missing visibilities are missing information, and deconvolution is therefore ill-posed and needs a constraint.
• Högbom's CLEAN (1974) supplies the constraint "the sky is a set of point sources" and proceeds greedily: find the peak, record a fraction (the loop gain) as a clean component, subtract a scaled-and-shifted dirty beam, repeat.
• We restore the component model by convolving it with a smooth clean beam (a Gaussian fit to the dirty beam's main lobe) and adding the residual back, giving a sidelobe-free image at the array's true resolution.
• We watched CLEAN pull our two point sources cleanly out of the sidelobes and reconstruct the extended blob, going from dirty image to a publishable restored image — the step that made aperture-synthesis imaging practical.

What's next¶

We now have the full classical synthesis-imaging pipeline in hand: visibilities → dirty image → CLEAN → restored image. The remaining chapters put real data through real tools. In Chapter 12 we will swap our NumPy toy for CASA's tclean and image an actual interferometric measurement set, seeing the major/minor-cycle machinery and the weighting and masking choices that turn Högbom's elegant core into a production imager.

---

📖 Glossary¶

Terms from this chapter — follow a link to look it up in the course glossary:

Aperture synthesis · Baseline · Beam · CASA · CLEAN · Deconvolution · Dirty beam / dirty image · Gain · Interferometer · Measurement Set · Point source · Sidelobe · Synthesised beam · Tapering / weighting