How much of the picture is the prior

A phase that survives what corrupts it ends by saying what a reconstruction does: the data determine a set of constraints, and the algorithm chooses among the pictures satisfying them by some criterion of smoothness or entropy or sparsity. That sentence is the whole of this essay, and it deserves a figure rather than a clause.

The underdetermination is not a matter of noise. An array samples the transform plane in arcs and leaves the rest of it untouched, and a sky brightness distribution has as many degrees of freedom as the field has resolution elements. Measure a few dozen of them and the rest are free — not poorly measured, free.

The family, stated exactly #

Let the sky be a vector $\mathbf{f}$ and the measurement a linear operator $A$ that takes it to a handful of visibilities $\mathbf{v} = A\mathbf{f}$. Any two skies differing by something in the null space of $A$ — anything $\mathbf{g}$ with $A\mathbf{g} = 0$ — produce identical data.

The set of skies consistent with a measurement is therefore

$$\{\mathbf{f}_0 + \mathbf{g} : A\mathbf{g} = 0\},$$

an affine subspace whose dimension is the number of unmeasured degrees of freedom. For the figure above that is 118 out of 128. For a real millimetre array imaging a compact source it is of the same order: a few hundred independent visibility measurements against an image with tens of thousands of pixels.

Nothing about this is specific to interferometry. It is the shape of every inverse problem with fewer measurements than unknowns, and the standard treatment is the same everywhere: add a penalty, minimise it over the family, and report the minimiser.

The choice of penalty is where the physics is supposed to enter, and where it usually does not.

Smoothness — minimise the total squared gradient — says the sky has no structure finer than it needs. It produces images with no sharp edges, which is a fair description of a nebula and a poor one of a set of point sources — the same choice an adaptive-optics pipeline makes when it decides how many components of a speckle field to remove.

Sparsity — minimise the total absolute brightness — says the sky is mostly empty. It produces images made of isolated compact components, which is a fair description of a set of point sources and a poor one of a nebula.

Maximum entropy — maximise $-\sum f_i \ln f_i$ — says the sky is as featureless as the data permit, with positivity built in. It sits between the two, and its output depends on the assumed default image in a way that is easy to overlook.

And the classical algorithm of radio astronomy is none of these formally, and is a sparsity prior in practice: it iteratively subtracts scaled copies of the point-source response at the brightest remaining position, so the model it builds is a sum of delta functions by construction.

Where the reconstructions agree, and why #

The important thing about the hero figure is not that the curves differ; it is where.

They agree in the measured frequencies, because they must. They agree at the location of the two compact sources, because a source’s position is carried by the phase of every measured visibility and is the best-determined quantity in the whole problem. And they agree, roughly, on the total flux, because the zero-frequency visibility was measured.

They disagree about whether the two compact sources are separated by a gap or joined by a bridge, about how much broad emission there is, and about the sharpness of everything — all of which are statements about the frequencies that were not sampled.

A reconstruction is trustworthy in proportion to how much of the claimed feature is in the measured frequencies, and the useful discipline is to say which those are. A feature at a spatial scale the array sampled well is a measurement; one at a scale it sampled once is a suggestion; one at a scale it never reached is the regulariser.

What the honest tests are #

The field’s practice has converged on three, and none of them is “the algorithms agree”.

Fit the visibilities, not the image. Where the quantity of interest is a number rather than a picture — the diameter of a ring, the separation of a binary, the flux ratio — it is fitted directly to the measured visibilities with a parameterised model, and no reconstruction enters. That is how the black hole shadow’s diameter was obtained: the published angular size is a fit to the data and the image is an illustration of it. The distinction matters because a fitted parameter has a likelihood and a reconstructed feature does not.

Run the pipeline on a known sky. Simulate an observation of a source whose truth is stipulated, with the same array, the same coverage and the same noise, reduce it exactly as the real data were reduced, and see what comes out. This is the same discipline that a high-contrast imaging pipeline applies by injecting fake planets, and it converts an unknowable systematic into a measured one.

And survey the priors deliberately. Rather than running several algorithms and hoping they agree, run one algorithm over a grid of its own parameters and over several priors, and report the range. The published black-hole images were produced this way: thousands of reconstructions across a parameter survey, with the ones that failed the synthetic-data test discarded, and the released image an average over the survivors.

The third of those is the one worth understanding, because it inverts the usual reading. Agreement between two algorithms is weak evidence — two sparsity-based methods will agree with one another and both be wrong about a smooth source. What the survey establishes is the range of images consistent with the data under any prior the synthetic tests do not rule out, and the claim is about that range.

The algorithm that looks like an image and is a model #

The oldest and still most-used reconstruction in radio astronomy deserves setting out, because its structure explains a habit that looks odd from outside.

Begin with the dirty image. Find its brightest pixel. Subtract a small fraction — a tenth, say — of the dirty beam, scaled and centred there, and record what was subtracted as a component. Repeat, thousands of times, until the residual looks like noise. The output is a list of components, not a picture.

To make a picture, the components are convolved with a clean, sidelobe-free beam of the same width as the dirty one’s main lobe, and the residual is added back. The published image is therefore a model deliberately blurred to the resolution the data support, plus whatever the algorithm could not account for.

Three features of that follow directly, and all three are unobvious.

The restoring beam is a choice. Making it narrower than the main lobe produces a sharper image with no more information in it, and the practice — which has a name, superresolution, and a bad reputation — is legitimate only when the model’s components are genuinely better localised than the beam, which is a question about the source rather than about the algorithm.

The residual is in different units. It is added back unconvolved, so a restored image is a sum of two things with different effective resolutions, and a faint extended feature in it may be residual rather than model.

And the components are not physical. A smooth source is represented as a dense forest of point components that happen to sum to something smooth, and reading individual components as sources is a well-documented way to invent structure. The extensions that handle extended emission replace the point component with a set of scaled shapes, which helps and introduces a shape as a prior.

The reason this algorithm has lasted fifty years despite not minimising anything in particular is that it is fast, it is stable, and its implicit prior — the sky is a set of compact things — matches the sources radio astronomers mostly look at. On the sources it does not match, it is known to fail in a characteristic way, producing a stripey residual at the scale of the missing short baselines.

Why the problem is unusually sharp here #

Every field with an inverse problem has this difficulty, and two things make it worse in interferometry than in most.

The data are extraordinarily few. A millimetre array of eight sites yields 28 baselines, of which perhaps 20 give useful signal-to-noise, tracked over a few hours. That is a few hundred independent complex numbers to constrain an image, against a medical scanner’s millions of projections or an optical telescope’s every pixel — and an optical interferometer has fewer still.

And the phases are mostly missing. Closure relations recover $(N-2)/N$ of the phase information, so a small array is working with a fraction of the constraint even on the frequencies it did sample. An eight-site array recovers three-quarters; a three-element optical interferometer recovers a third.

Put together, a millimetre-VLBI reconstruction has fewer constraints per degree of freedom than almost any imaging problem in science, and correspondingly more of the picture is supplied by the choice of penalty. That is not a criticism of the images; it is the reason the surveys are run and the reason the published uncertainties on image-derived quantities are so much larger than the uncertainties on fitted ones.

The missing spacings #

There is one hole in the coverage that no amount of observing can fill, and it produces the most common quantitative error in interferometric imaging.

Two antennas cannot be closer together than their own diameters. So the shortest baseline an array can form is of order one dish width, and every spatial frequency below that — every structure broader than the corresponding angular scale — is unmeasured. An interferometer is not merely insensitive to smooth large-scale emission; it is exactly blind to it, and its zero-spacing visibility, which is the source’s total flux, is one of the frequencies it cannot reach.

The consequence is a systematic with a name and a sign. An extended source observed by an array alone comes out with a negative bowl around it: the reconstruction, having no information about the broad component, supplies a smooth negative background to make the measured visibilities work out. The source’s total flux is then underestimated, sometimes by a large factor, and the shape of what is left is distorted at the largest scales.

The fix is to measure the missing frequencies with a different instrument — a single large dish, which samples exactly the short spacings an array cannot — and combine the two. That is standard practice and it is not free: the two data sets have different calibrations, different resolutions and different noise, and how to weight them together is another choice made outside the data.

The largest scale an array can see is set by its shortest baseline and the smallest by its longest, and a source larger than the first is measured with a hole in the middle of its transform. That hole is in the null space like any other, with the difference that it is at the frequencies carrying most of the flux.

What was actually measured #

Complex visibilities: for each pair of sites and each integration, a correlation amplitude and a phase, with the phase usually unusable on its own and entering through closure quantities.

Three things about that data set are worth naming because they enter every reconstruction and none of them is in the figures.

The amplitudes need a calibration the phases do not. A visibility amplitude is a flux, and a flux requires knowing each antenna’s gain — which drifts with elevation, weather and receiver temperature. Closure amplitudes remove it at the cost of information, exactly as closure phases do.

The noise is not Gaussian in the quantity being fitted. A visibility’s real and imaginary parts have Gaussian errors; its amplitude does not, because an amplitude is a positive quantity formed from two noisy numbers and is biased upward at low signal-to-noise. Fitting amplitudes without correcting for that bias systematically inflates a source’s compact flux.

And the coverage is not a property of the array. It depends on the source’s declination, on how long it was tracked, and on which sites had weather. Two observations of the same source with the same array can have substantially different null spaces, and therefore different sensitivity to the prior.

Where the model stops #

The figures are one-dimensional. A real image is two-dimensional, the null space is correspondingly larger, and the penalties behave differently — a two-dimensional smoothness prior has a direction as well as a strength, and an elongated beam imposes structure that the reconstruction has to be prevented from believing.

Positivity is not imposed here and is imposed in practice. The sky’s brightness cannot be negative, and that constraint is genuinely informative: it removes a great deal of the null space, which is why the maximum-entropy methods that build it in outperform the unconstrained ones. The reconstructions drawn here do not use it, which makes their disagreement larger than a real pipeline’s would be — the figure overstates the size of the effect and not its existence.

The truth is stipulated. Everything here is a synthetic sky whose answer is known, which is what allows the three reconstructions to be judged. In a real observation there is no such curve to draw, and the comparison that can be made is between reconstructions rather than against the sky.

And the noise is absent. With noise, “reproducing every measured visibility exactly” is the wrong objective — an image that fits noisy data exactly has fitted the noise. Real reconstructions minimise a penalty subject to a statistical fit, and how tightly to fit is another parameter with a prior in it.

The generalisation #

The structure is the one every underdetermined inverse problem has, and the useful form to carry is about which questions a data set can answer rather than how well it answers them.

A measurement operator has a null space. Anything in it is invisible, exactly and forever, and no improvement in the instrument’s sensitivity touches it. So the first question about any reconstruction is not how accurate but what is in the null space — and the second is whether the quantity being claimed has a component in it.

A number fitted to the data has an error bar; a feature read off a reconstruction has a prior. Those are different kinds of statement and the difference is not one of degree. The discipline that follows is to reduce every claim to a parameter of a model wherever possible, and to treat the picture as a way of seeing what to fit rather than as the result.

There is a second reading that runs the other way and is worth keeping beside the first. The prior is not a contaminant to be removed; it is what makes the problem soluble at all, and it encodes real knowledge — the sky is non-negative, it is not infinitely spiky, it does not have structure finer than the physics allows. An inverse problem with no prior has no answer. The objection is never to having a prior but to not saying what it was, and the fields that handle underdetermination well are the ones that publish the penalty alongside the picture.

Still open: the position that was thrown away #

What comes next takes the other branch. Closure quantities are immune to the atmosphere because they are blind to position, and the two are the same degree of freedom — so an image made from them is correct in structure and floating on the sky. Recovering the position means giving the immunity up deliberately, referencing the phase to a nearby calibrator and measuring the atmosphere rather than cancelling it.

Beside it lies the question treated as settled here: what the array should be, given that coverage rather than algorithm is the cure. The design of an array is the choice of a null space, and it is made decades before anybody knows what will be observed with it.