Local cosmography and the Wiener filter

The local cosmic neighbourhood has recently been mapped in spectacular fashion by the Cosmicflows research programme, yielding papers in Nature, and an article in Scientific American. The images generated are stunning.

Perspective view of the X-Y equatorial plane of our cosmic neighbourhood in supergalactic coordinates. The density of matter is represented by colour contours, deep blue regions indicating voids, red indicating zones of high density. The velocity streams are represented in white, with individual galaxies as spheres.

There's also an interesting mathematical back-story here because the work has been underpinned by techniques developed over 20 years ago by Yehuda Hoffman and colleagues. Drawing upon the exposition provided by Hoffman, let's take a look at these methods, beginning with Gaussian fields.

Gaussian fields

A random field is a field in which the value at each point is sampled from a univariate probability distribution, and the values at multiple points are sampled from a multivariate distribution. Typically, the sampled values at different points are spatially correlated, with a degree of correlation that declines with distance.

A Gaussian field is a type of random field in which the distribution at each point is Gaussian, and relationship between the values at different points is given by a multivariate Gaussian distribution. The properties of a Gaussian field are specified by its covariance matrix. Assuming a Gaussian field $\mathbf{s}$ of zero mean, this is denoted as:
$$
\mathbf{S} = \langle \mathbf{s} \mathbf{s}^T\rangle
$$ Treating $\mathbf{s}$ as a vector of random variables, this expression is understood as the 'outer-product' of the column vector $\mathbf{s}$ with its transpose row vector $\mathbf{s}^T$:
$$
\mathbf{S} = {\langle s_i s_j \rangle}
$$Convolution

Given the Gaussian field $\mathbf{s}$, the generation of the measured data-set $\mathbf{d}$  is represented as follows:
$$
\mathbf{d} = \mathbf{Rs}  + \mathbf{\epsilon}
$$ $\mathbf{R}$ represents the characteristics of the measuring instrument, and $\epsilon$ represents measurement noise. The covariance matrix of the noise is denoted as follows:
$$
 \mathbf{N} = \langle \mathbf{\epsilon} \mathbf{\epsilon}^T\rangle
$$ The transformation $\mathbf{R}$ is typically a convolution. In astronomical measurements of the sky it is often referred to as the Point Spread Function. It specifies how the energy detected to belong to a particular pixel of the sky is actually a weighted sum of the energy belonging to a range of adjacent pixels. Similarly, in spectroscopy it specifies how the energy detected at each particular wavelength is actually a weighted sum of the energy belonging to adjacent wavelengths, thereby smearing out the energy spectrum.

Given the true spectrum $f(\lambda)$, and the convolution kernel $g(\lambda - \lambda_1)$, the measured spectrum $I(\lambda)$ is defined by the convolution:
$$
I(\lambda) = f * g = \int f(\lambda_1) g(\lambda - \lambda_1) d\lambda_1
$$

Wiener filter

Stochastic filtering provides a collection of techniques for constructing estimates of the true values of a quantity or field from sparse and noisy measurement data. The Wiener filter is one such technique. It provides a transformation $\mathbf{F}$ that maps a measured dataset $\mathbf{d}$ into an estimate of the true field:
$$
\mathbf{s}^{WF} = \mathbf{F} \mathbf{d}
$$ The discrepancy between the true field and the estimated field is called the residual:
$$
\mathbf{r}  = \mathbf{s} - \mathbf{s}^{WF}
$$ The residual possesses its own covariance matrix $\langle \mathbf{r}\mathbf{r}^T \rangle$. The Wiener filter is defined so that it minimizes the covariance of the residual.

The Wiener filter possesses another property which makes it the natural choice of filter for cosmography: Given the prior distribution of the true field $\mathbf{s}$; given the relationship between the true field and the measured values; and given the measured values, a posterior Bayesian distribution $p(\mathbf{s}|\mathbf{d})$ can be obtained over the true field. In the case where the true field is a Gaussian field, and the noise is also Gaussian, the Wiener filter picks out the mean value of the posterior distribution.

The Wiener filter is given by the expression:
$$
\mathbf{F} = \mathbf{SR}^T (\mathbf{RSR}^T + \mathbf{N})^{-1}
$$

Constrained realizations

The method of constrained realizations goes back to a paper by Yehuda Hoffman and Erez Ribak in 1991 ('Constrained realizations of Gaussian fields - A simple algorithm', Ap. J. Lett., vol 380, pL5-L8). It's based upon a simple but crucial point: the application of the Wiener filter to a measured dataset will produce an estimated field which includes the covariance of the residual: $$\mathbf{s}^{WF} = \mathbf{s} - \mathbf{r}$$

Hence, the estimated field will be smoother than the true field. The idea proposed by Hoffman and Ribak is simple: to generate a realistic realization of the true field, you need to add a sampled realization from the residual field.

Their method works as follows:

(i) Generate a random realization $\tilde{\mathbf{s}}$ of the true field. (The tilde here indicates a realization).

(ii) Generate a realization of the measurement noise, and apply the measurement transformation to $\tilde{\mathbf{s}}$ to yield a random dataset realization: $$\mathbf{\tilde{d}} = \mathbf{R\tilde{s}}  + \mathbf{\tilde{\epsilon}}$$

(iii) Apply the Wiener filter to $\mathbf{\tilde{d}}$ to create an estimate of the true field realization: $\mathbf{F} \mathbf{\tilde{d}} $.

(iv) Generate a realization of the residual:

$$\mathbf{\tilde{r}} = \tilde{\mathbf{s}} - \mathbf{F} \mathbf{\tilde{d}} $$

(v) Add the realization of the residual to the estimated field which has been obtained by applying the Wiener filter to the actual measured dataset: 

$$\eqalign{\mathbf{s}^{CR} &=  \mathbf{\tilde{r}} + \mathbf{Fd} \cr &=\tilde{\mathbf{s}}+ \mathbf{F} (\mathbf{d-\tilde{d}})}$$

Reconstruction of continuous fields from finite data

Here the Gaussian field $\mathbf{s}$ is assumed to be a function $f(\mathbf{r})$, so that the covariance matrix becomes the auto-correlation function $\xi$:

$$
\mathbf{S} = \langle f(\mathbf{r}_i) f(\mathbf{r}_j)\rangle =  \xi(\mathbf{r}_i,\mathbf{r}_j)
$$ Assume that measurements are made of this field at a finite number of points, ${F_i}$. If there is no convolution operation, so that $\mathbf{R} = \mathbf{I}$, then in this case the combined application of the Wiener filter and constrained realization yields the following expression:
$$
 f(\mathbf{r})^{CR} = \tilde{f}(\mathbf{r})+ \xi(\mathbf{r}_i,\mathbf{r}) ( \xi(\mathbf{r}_i,\mathbf{r}_j) + \mathbf{N}_{ij} )^{-1}(F_j-\tilde{f}(\mathbf{r}_j))
$$ Repeated indices are summed over here.

Cosmicflows

The Cosmicflows research programme applied these techniques to a finite collection of sparse and noisy galaxy recession velocities to reconstruct the full 3-dimensional peculiar velocity field in our local cosmic neighbourhood. (The latter is defined to be a length-scale of the order 100 Mpc). 'Peculiar' velocities in this context are those which remain after the recession velocity due to cosmic expansion has been subtracted.

Our cosmic neighbourhood,within a cube 200Mpc on each side. The Milky Way is located at the origin of coordinates. There are three isodensity colour contours. The velocity flows are depicted as black lines.

Assuming that cosmic structure formation has been seeded by perturbations sampled from a Gaussian random field, and assuming that the expansion of the universe followed the Lambda-Cold Dark Matter ($\Lambda$CDM) model, the combination of the Wiener filter and Constrained Realization was applied to a sample of galaxy recession velocities along the line-of-sight. The density field in our neighbourhood was then inferred from the reconstructed velocity field.

Some of the most striking images produced by the Cosmicflows team are those which depict the streamlines of the peculiar velocity field in our cosmic neighbourhood. These reveal regions in which the streamlines converge, due to concentrations in the density of matter such as the Great Attractor, and regions in which the streamlines diverge from cosmic voids. In particular, the Cosmicflows team have attempted to identify the surface of divergent points that surrounds us. They define the volume within that watershed as our local galaxy Supercluster, which they name 'Laniakea'.

The Cosmicflows programme assumes only linear deviations from homogeneity and isotropy. This assumption preserves the Gaussian nature of the peculiar velocity field and the density perturbation field. But it does have the consequence that the velocity field is irrotational. i.e, it possesess zero vorticity. Hence, the images of cosmic streamlines contain watersheds, but no whirlpools. Whilst the linear theory should prevail on large scale, the non-linearity should dominate on smaller scales. The Cosmicflows team claim that the non-linearity should only dominate on the Mpc length-scale, and that their approach is therefore valid.

A slice through the equatorial plane in the supergalactic coordinate system. The boundary of the Laniakea supercluster is depicted as an orange line. The diagram is shaded with density contours, deep blue regions indicating voids, and red regions indicating zones of high density. The velocity streams within the Laniakea basis of attraction are represented in white.