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.
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.

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}

$$

__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

$$

$$

\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.

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.