The OED formalism — Fisher information, the additive backbone, and sky projection

The OED formalism — the machinery every example reuses¶

Every worked example in this section — the anisotropy design, the dynamical-mass depth knob — runs on the same engine. This page builds that engine once, from first principles, so the example pages can spend their words on the science instead of re-deriving the algebra. If you have just landed here from one of them, this is the page that explains why a pre-data error bar is even computable.

Why design an observation at all?¶

Telescope time is the scarcest resource in observational astronomy. A night on an 8-metre is worth tens of thousands of dollars; a Gaia-quality proper motion or a high-resolution spectrum of a faint cluster star is expensive in a way that the analysis afterwards is not. And yet the decision of what to observe — how many stars, how far out, how faint, in which channel — is usually made by intuition and heritage: “measure the brightest stars near the centre, that’s where the signal is.”

Sometimes intuition is right. Often it is exactly backwards. The information about a parameter is not where the signal is largest; it is where the signal is most sensitive to that parameter. Those can be very different places, and the gap between them is telescope time thrown away.

Optimal experimental design (OED) replaces the guess with a computation. The key realisation is that the precision you will achieve is a known function of the observing strategy — known before you collect a single photon — through the Fisher information. So you can treat the observing strategy as a free variable and maximise the information you expect to extract. The output is not just “a good plan”; it is a quantified plan: “this allocation reaches the same precision with $3.66\times$ fewer stars,” or “going one magnitude deeper buys you nothing.” OED is rare in astronomy because almost no simulator can compute $\partial(\text{information})/\partial(\text{observing strategy})$ — and that is exactly what a differentiable forward model makes possible.

The foundation: Fisher information and the Cramér–Rao bound¶

Suppose your data $\mathbf{d}$ depend on parameters $\boldsymbol\theta$ through a likelihood $\mathcal{L}(\mathbf{d}\mid\boldsymbol\theta)$, and write the log-likelihood $\ell = \ln\mathcal{L}$. Near the truth, $\ell$ is a hill in parameter space, and its curvature tells you how sharply the data pin the parameters: a sharp peak means a tight measurement, a flat ridge means a degeneracy. The Fisher information matrix is exactly that curvature, averaged over data realisations:

Its meaning is delivered by the Cramér–Rao bound: no unbiased estimator can do better than

So $F^{-1}$ is the best achievable covariance, and $\sqrt{(F^{-1})_{ii}}$ is the tightest $1\sigma$ error bar you can hope to put on $\theta_i$ — computable from the model alone, before observing. That is the hinge of the whole method: the error bar is a function of the design, so you can optimise it.

For data with Gaussian errors — our binned dispersions, with measurement uncertainty $\delta\sigma$ — the Fisher matrix has a simple, additive form. If $g(\boldsymbol\theta)$ is the model prediction for a datum,

Two features of Equation do all the work below. First, it is a sum over independent data — Fisher information adds. Second, each term is the sensitivity $\partial g/\partial\theta$ squared, divided by the noise. A datum is informative not when $g$ is large, but when $g$ moves as $\theta$ moves.

The dimensionless $\ln\theta$ metric¶

Cluster parameters span enormous dynamic ranges: a total mass $M\approx10^5\,\Msun$ sits next to a half-mass radius $r_h\approx3$ pc in the same vector. Built in raw units, the Fisher is wildly ill-conditioned ($\mathrm{cond}\approx1.7\times10^9$) and the A/D criteria below are not scale-invariant — $\mathrm{tr}\,F^{-1}$ would be dominated by whichever parameter happens to carry the largest unit magnitude, not the largest fractional uncertainty.

The additive, design-linear backbone¶

This is the structural trick that makes the whole optimisation cheap. Start from the per-datum noise: a velocity dispersion measured from $n$ stars, each with intrinsic dispersion $\sigma$ and measurement error $\epsilon$, has a sampling uncertainty

The factor of $2n$ (not $n$) is the standard error on a dispersion, not a mean — a dispersion is a second moment, so it converges more slowly. More stars in a bin ($n\uparrow$) sharpen that bin’s datum; that is the design’s only lever on the noise.

The design and the budget. A design allocates a fixed budget $N_{\rm total}$ across the cells (projected-radius bin $b$ × kinematic channel $c$). We parametrise it on the softmax simplex so it stays a valid, positive allocation under unconstrained optimisation,

with the logits $\mathbf z$ as the design variables. Substituting Equation into Equation and grouping by cell, the design Fisher becomes a budget-weighted sum of design-independent per-star blocks:

Each $M_{b,c}$ is a $3\times3$ rank-1 matrix — the information one star in that cell-and-channel carries about $\boldsymbol\theta$ — and it does not depend on the design. So the Jacobian $J_{b,c}$ is computed once, by a single reverse-mode jacrev through the differentiable forward model at the truth. After that, the design enters Equation only through the scalar weights $n_{b,c}$:

Optimising the observing strategy is $F=\sum n\,M$ → invert a $3\times3$ → read one element. The expensive projection is never re-differentiated inside the optimiser loop, and the gradients $\partial(\text{criterion})/\partial(\text{design})$ are pure linear algebra over the precomputed blocks.

We use jacrev (reverse mode) because it is the supported, tested AD path for all profiles and stays correct if a King/Michie mock is ever swapped in (those equilibrium-solver profiles hit a custom_vjp ODE with no forward-mode rule). On the analytic-density Plummer path used in the examples there is no ODE, so forward-mode would also work — jacrev is the robust choice, not a forced one.

What “optimal” means: the c / D / A criteria¶

The same $F=\sum n\,M$ supports several notions of “best.” The choice of criterion is the choice of what you are optimising for:

c-optimality is the headline of both examples: it minimises the variance of one target parameter $\theta_\star$ after marginalising over the nuisances — exactly the $(\theta_\star,\theta_\star)$ element of the full inverse, which already accounts for how the target trades off against the others. D-optimality maximises the determinant of $F$, shrinking the volume of the joint confidence ellipsoid over all parameters. A-optimality minimises the trace of $F^{-1}$, the average variance. They are genuinely different objectives, and under parameter degeneracies they put stars at different radii.

The projection geometry: $\sigma_r(r)$ into the three sky observables (B&M82)¶

Both examples observe the cluster through the same three kinematic channels, and the geometry that connects the cluster’s internal velocity ellipsoid to what a telescope records is worth getting exactly right — it trips up nearly everyone the first time. In a spherical system each star’s velocity splits into a radial part (toward or away from the cluster centre) and two tangential parts, with dispersions $\sigma_r$ and $\sigma_t$. The Binney anisotropy parameter measures their asymmetry,

which is 0 for isotropic orbits and rises toward 1 for radially-biased ones.

The word “radial” means three different things here, and conflating them is the classic mistake:

The radial velocity you measure spectroscopically is radial to the observer — the line-of-sight Doppler shift — not the 3-D motion toward the cluster centre. A star at projected radius $R$ and line-of-sight depth $z$ sits at 3-D radius $r=\sqrt{R^2+z^2}$, and its local radial direction $\hat r$ makes an angle $\psi$ with the sight line, $\cos\psi = z/r$. The observed line-of-sight velocity dispersion at that point is therefore a projection that mixes the radial and tangential dispersions:

using $\sigma_t^2=(1-\beta)\sigma_r^2$ and $z^2+R^2=r^2$. Only for a sight line through the centre ($R\to0$, $\psi\to0$) does $\sigma_{\rm los}=\sigma_r$; everywhere else the LOS velocity is a $\beta$-weighted blend. Integrating Equation along the sight line, weighted by the density $\rho(r)$, gives the observed profile. Doing the same for the two plane-of-sky (proper-motion) components yields the full Binney & Mamon (1982) projection — three distinct observables of the same underlying $\sigma_r(r)$ and $\beta(r)$:

The $\sigma_r^2(r)$ inside these integrals is itself the solution of the anisotropic Jeans equation; for the Osipkov–Merritt law $\beta_{\rm OM}(r)=r^2/(r^2+r_a^2)$ Merritt, 1985 it has the closed integrating-factor form

evaluated by the packaged differentiable project_dispersion. The point of all this machinery is one sentence:

In the isotropic limit ($\beta=0$) the three kernels in Equation collapse to 1 and $\sigma_{\rm los}=\sigma_{{\rm pm},R}=\sigma_{{\rm pm},T}$ — the anisotropy lives only in the ratios between channels.

The tangential proper motion is the cleanest probe: $\sigma_{{\rm pm},T}^2\propto (1-\beta)$, so as orbits turn radial ($\beta\to1$) the tangential PM collapses while the LOS and radial-PM channels stay up. And because $\beta(r)$ grows outward, the channels disagree most in the outskirts. A single $\sigma_{\rm los}(R)$ profile cannot break the tie — a cold outer LOS profile can mean “low $\sigma_r$” or “high $\beta$.” Proper motions break it. Notice too that the total mass $M$ enters every one of these channels through $GM(<s)$ in Equation: $\sigma_r^2\propto GM$, so a dispersion measured anywhere weighs the cluster — the lever the dynamical-mass example pulls.

References¶

The Fisher-information / Cramér–Rao foundation is standard estimation theory (see e.g. Tegmark, Taylor & Heavens 1997 for the astrophysical forecasting form). The Osipkov–Merritt anisotropy law $\beta(r)=r^2/(r^2+r_a^2)$ is Merritt (1985); the line-of-sight projection of an anisotropic spherical model into $\sigma_{\rm los}$, $\sigma_{{\rm pm},R}$, $\sigma_{{\rm pm},T}$ is Binney & Mamon (1982, MNRAS 200, 361). The differentiable project_dispersion forward model is documented on the velocity-DF kinematics pages.