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.
- Implement Högbom's CLEAN yourself, in flat NumPy, on a real dirty image and
watch the sources emerge from the sidelobes -- then promote the loop to the
packaged, tested
jansky.interferometry.hogbom_cleanfor heavier runs. - See, qualitatively, how CASA's
tcleangeneralises 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:
- Find the peak. Locate the pixel $(i, j)$ of largest absolute value in the current residual image (which starts equal to the dirty image).
- 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)$.
- 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.
- 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.
You are about to write exactly these four steps yourself, in the code below.
Later in the chapter you will meet jansky.interferometry.hogbom_clean -- the
identical loop packaged in src/jansky/interferometry.py -- for the heavier runs.
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.
Back of the envelope¶
Before any code: how many CLEAN iterations does the algorithm actually need? Suppose the dirty image below is dominated by one bright point source at a normalised peak of 1 Jy/beam, and we run Högbom's CLEAN with the standard loop gain $g = 0.1$. Every iteration removes a fraction $g$ of whatever peak remains, so after $n$ iterations the leftover peak is a fraction $(1-g)^n$ of the original.
The question: roughly how many iterations does it take to clean that source down to 0.1% of its peak flux -- a stand-in for a typical few-$\sigma$ noise floor?
Rules of the napkin: decades matter, factors of two don't. Work out the algebra on paper first -- solve $(1-g)^n = \text{target}$ for $n$ -- then write it as one line of code below.
import math
S_peak = 1.0 # Jy/beam -- normalised peak of the dominant source in the dirty image
gain = 0.1 # loop gain: CLEAN removes this fraction of the current peak each iteration
target_frac = 1.0e-3 # residual target, as a fraction of the peak (a few-sigma noise floor)
# Your one line: solve (1 - gain)**n = target_frac for n. (hint: math.log(x, base))
n_iters_guess = None # <- replace None with an expression in gain and target_frac
To grade an envelope estimate we only care about the exponent -- the course's
one-line checker, jansky.envelope.check (introduced in Chapter 1), compares your
guess to the worked answer by decades rather than decimal places.
from jansky.envelope import check
# The answer is worked in the reveal below -- commit to your own number first.
# The check stores only its base-10 logarithm, so a stray glance can't spoil it.
check(n_iters_guess, expected_log10=1.8195, name="iterations to 0.1% residual", units="iterations")
[iterations to 0.1% residual] no guess yet — fill in the cell above, then re-run this one.
False
Reveal — the envelope answer
Each iteration multiplies the remaining peak by $(1-g)$, so after $n$ iterations only a fraction $(1-g)^n$ is left. Solve for $n$:
$$ (1-g)^n = \text{target} \quad\Longrightarrow\quad n = \frac{\ln(\text{target})}{\ln(1-g)} = \frac{\ln(10^{-3})}{\ln(0.9)} \;\approx\; 66. $$
n_iters_guess = math.log(target_frac, 1 - gain) # = 65.6 -> ~66 iterations
Sixty-six iterations, for a residual a thousand times fainter than the peak. Hold that number: the actual run below (a real, multi-source dirty image, not an idealised single point) crosses this same 0.1%-of-peak threshold within well under a factor of 2 of this estimate -- even though the napkin math ignored everything except the brightest source.
Where the envelope leaks (and where this chapter patches it):
- We assumed one dominant point source taking the peak every single iteration. Real dirty images have several sources and sidelobe structure competing for the peak -- see the real (multi-source) run in "Running CLEAN" below.
- We tracked only the peak pixel value, not whether CLEAN puts the flux in the right place. "Did CLEAN actually recover the sources?" checks the positions and relative fluxes, not just how fast the peak falls.
- We picked an arbitrary target (0.1% of peak) as a stand-in for "the noise floor." Exercise 2 shows what happens if you keep iterating well past that point -- you start "cleaning the noise."
Turn the knob¶
Exercise 1 later in this chapter sweeps the loop gain from 0.05 up to 0.8. Same source, same 0.1%-of-peak target -- but crank the gain to $g = 0.5$, five times more aggressive. Reuse your expression.
gain_aggressive = 0.5 # a much more aggressive loop gain
n_iters_aggressive_guess = None # <- your expression again, with the new gain
check(
n_iters_aggressive_guess,
expected=10.0,
name="g=0.5, iterations to 0.1% residual",
units="iterations",
)
[g=0.5, iterations to 0.1% residual] no guess yet — fill in the cell above, then re-run this one.
False
Reveal — the scaling answer
$$ n = \frac{\ln(10^{-3})}{\ln(1 - 0.5)} \;\approx\; 9.97 \;\approx\; 10. $$
n_iters_aggressive_guess = math.log(target_frac, 1 - gain_aggressive) # = 9.97 -> ~10
Five times the loop gain buys roughly 6-7 times fewer iterations to the same residual target -- the relationship is logarithmic, not linear, so a bigger gain has real but diminishing returns. The catch, which Exercise 1 explores directly, is that a large gain claims flux faster but is more likely to over-fit noise and sidelobe structure instead of genuine sky structure -- the speed-vs-stability trade CLEAN always makes.
Setup¶
We reuse the same helpers as Chapter 8 — jansky.interferometry for the
synthesis machinery (uv-coverage, dirty beam, dirty image) and
jansky.plotting for consistent figures. CLEAN itself we build by hand below,
one step at a time -- the packaged interferometry.hogbom_clean only reappears
once we've written that loop ourselves, for the heavier exploratory runs later
in the chapter.
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
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
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.
# 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))
sky peak = 1.0 beam peak = 1.0 (normalised) dirty peak = 0.087
# 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 now do exactly what Högbom described, step by step, in the open: find the peak of the residual, record a fraction of it (the loop gain) as a clean component, subtract a beam scaled and shifted onto that pixel, and repeat. This is the whole algorithm -- about a dozen lines of NumPy -- and it stays inline here because deconvolution is this chapter's subject. We keep three running arrays as we go:
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.
We use a conservative loop gain of $g = 0.1$ and a generous budget of 500 iterations.
gain = 0.1 # loop gain: subtract 10% of each peak per iteration
n_iter = 500 # maximum number of CLEAN iterations
# 1. Start the residual as the dirty image, with an empty component model.
residual = dirty.copy()
model = np.zeros_like(dirty)
components = [] # (row, col, flux), in subtraction order
beam_center = np.array(beam.shape) // 2 # the dirty beam is centred on the grid's middle pixel
for _ in range(n_iter):
# 2. Find the peak of the current residual -- CLEAN's guess at "a real source is here".
idx = np.unravel_index(np.argmax(np.abs(residual)), residual.shape)
peak = residual[idx]
# 3. Record a fraction (the loop gain) of that peak as a clean component.
flux = gain * peak
model[idx] += flux
components.append((int(idx[0]), int(idx[1]), float(flux)))
# 4. Subtract a copy of the dirty beam, scaled by flux and shifted onto the peak pixel --
# this removes the component AND all of the sidelobes it would have thrown.
shift = (idx[0] - beam_center[0], idx[1] - beam_center[1])
shifted_beam = np.roll(np.roll(beam, shift[0], axis=0), shift[1], axis=1)
residual = residual - flux * shifted_beam
# 5. Repeat -- this loop always runs to n_iter (no stopping threshold set here).
print(f"CLEAN found {len(components)} components in {n_iter} iterations")
print(
f"residual peak = {np.abs(residual).max():.4f} "
f"(dirty peak was {np.abs(dirty).max():.4f})"
)
print("brightest few components (row, col, flux):")
for row, col, flux in components[:5]:
print(f" ({row:3d}, {col:3d}) flux = {flux:.4f}")
CLEAN found 500 components in 500 iterations residual peak = 0.0001 (dirty peak was 0.0865) brightest few components (row, col, flux): (120, 110) flux = 0.0087 (120, 110) flux = 0.0078 (120, 110) flux = 0.0070 (120, 110) flux = 0.0063 (120, 110) flux = 0.0057
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.
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(model, sigma=sigma_clean)
restored = restored_model + residual
print("restored image peak =", round(float(restored.max()), 3))
clean-beam sigma ~ 1.70 px (FWHM ~ 4.0 px) restored image peak = 0.005
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).
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 components]
comp_cols = [c for r, c, f in components]
show_image(
model, ax=axes[0, 1], title=f"2. Clean components ({len(components)} found)"
)
axes[0, 1].scatter(
comp_cols, comp_rows, s=10, facecolors="none", edgecolors="cyan", linewidths=0.6, alpha=0.5
)
show_image(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.
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 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(components, key=lambda t: t[2])
print(f"\nBrightest component is at ({peak_r}, {peak_c}); true bright source is at (120, 110).")
Recovered clean-component flux near each input source: bright (1.0 Jy) at (120, 110): summed component flux = 0.087 faint (0.6 Jy) at (150, 165): summed component flux = 0.052 Brightest component is at (120, 110); 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.
From napkin to package¶
Every step above -- find the peak, subtract the scaled beam, record the component --
was typed out by hand, because deconvolution is this chapter's subject. But every
later chapter that touches synthesis imaging (and the real CASA run in Chapter 12)
needs to run CLEAN on demand, at whatever gain and iteration budget the moment calls
for, without retyping the loop each time. So the course keeps one tested copy,
jansky.interferometry.hogbom_clean -- literally the loop you just wrote, wrapped in a
function with a docstring and a unit test. One assert proves it reproduces your model,
residual and component list exactly, and from here on we lean on it for the heavier
exploratory runs below.
from jansky.interferometry import hogbom_clean, CleanResult
packaged = hogbom_clean(dirty, beam, gain=gain, n_iter=n_iter)
assert np.allclose(packaged.model, model)
assert np.allclose(packaged.residual, residual)
assert packaged.components == components
print(f"inline CLEAN loop == jansky.interferometry.hogbom_clean -- promoted ({CleanResult.__name__}).")
inline CLEAN loop == jansky.interferometry.hogbom_clean -- promoted (CleanResult).
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.
# 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}")
gain residual peak n components 0.05 0.0001 500 0.10 0.0001 500 0.30 0.0000 500
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.
# 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,
tcleanalternates 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.