INTELLIGENT DESIGN: ON THE EMULATION OF COSMOLOGICAL SIMULATIONS

The very low expense of CMB temperature power spectrum calculations makes them convenient to use for testing and development of new emulation techniques. For this reason, we pursue emulation of a CMB temperature power spectrum calculator for our case study.

We set ourselves the goal of building an emulator with sufficient accuracy for analyzing the temperature power spectrum constraints that are expected from Planck,⁴ and to do so in the simplest possible manner with the fewest possible "simulations." To inform the emulator construction, we assume prior constraints on cosmological parameters of similar quality to those from WMAP5 (Nolta et al. 2009). We consider the following six cosmological parameters:

$$\begin{equation} {\theta } = (n_s, 100\, \theta ^{*}, \ln (10^{10} A), \omega _c, \omega _b, \tau)^{T}, \end{equation} \tag{ 1 }$$

where n_s is the scalar spectral index, θ* is the angular size of the sound horizon at last scattering (in radians), A is the amplitude of the primordial power spectrum defined at the pivot scale k = 0.002 Mpc⁻¹, ω_c ≡ Ω_ch² is the dark matter density, ω_b ≡ Ω_bh² is the baryon density, and τ is the optical depth. Our fiducial values for these parameters are given in the first column of Table 1 and are derived from a Markov Chain Monte Carlo (MCMC) chain using the WMAP5 likelihood.⁵ We use a Fisher matrix to compute constraints on these six cosmological parameters given a signal power spectrum computed with CAMB and a model for the noise power spectrum,

$$\begin{equation} C_{\ell }^{\rm noise} = \frac{e^{(\ell \,\sigma _{\rm beam})^2}}{w_{\rm noise}} \end{equation} \tag{ 2 }$$

with noise weights w_noise = 1.43 × 10¹⁴ for WMAP5, w_noise = 4.7 × 10¹⁶ for Planck, beam smearing σ_beam≡ FWHM/2.355, and FWHM = 0.00378 for WMAP5 (at 90 GHz) and 0.0021 for Planck (at 143 GHz).

We further include in our Fisher matrix a prior on a combination of optical depth and primordial power spectrum amplitude similar to what is provided by WMAP5 polarization data. Specifically, assuming the quantity Ae^−2τ has negligible error (since this combination is determined very well from the temperature power spectrum), we impose a Gaussian prior with width

$$\begin{equation} \frac{\sigma (A \tau ^2) }{A \tau ^2} \approx \frac{\sigma (e^{2\tau } \tau ^2) }{e^{2\tau }\tau ^2} \approx \frac{2(\tau +1)\sigma (\tau)}{\tau }, \end{equation} \tag{ 3 }$$

with σ(τ) = 0.017 from WMAP5,⁶ which comes from the combination of polarization data with temperature data and is the main contribution of the polarization data to parameter constraints. In Figure 1, we show the size of the power spectrum error bars for our WMAP5-like observations and those for our Planck-like observations. In Figure 2, we show the corresponding two-dimensional Fisher matrix contours.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We were sufficiently satisfied with the accuracy and simplicity of the emulator we developed, that we are releasing a version of it as a publicly available code, called Emu CMB, described in Appendix A. There is another fast and accurate code available (PICO; Fendt & Wandelt 2007a, 2007b) that has some advantages. Most notably, PICO is a more complete emulator that provides predictions for polarization as well as temperature power spectra and includes models with non-zero spatial curvature. Emu CMB's advantages are (1) its significantly simpler algorithm that has made it very easy for us to make available implementations in multiple languages with no need for external libraries, (2) the simple training algorithm can be reimplemented quickly (using our public code if desired) for revised CMB models, and (3) it extends to higher maximum multipole moment (ℓ_max = 5000). The latter quality is especially important for analysis of high-resolution data, such as from the Atacama Cosmology Telescope or the South Pole Telescope (SPT).

The design runs for this paper were all calculated using CAMB (2008 June release; Lewis et al. 2000) with accuracy settings accuracy_boost = 4, l_accuracy_boost = 4, l_sample_boost = 8, and without including gravitational lensing.⁷ With this accuracy setting, numerical integration errors are negligibly small. We modified CAMB to extract the power spectra only at multipoles where they are actually computed, skipping the default interpolation to all multipoles. For a maximum multipole of ℓ = 3000, this gives 448 power spectrum values using the default multipole spacing in CAMB.

Here we review two design procedures: that used for CC and that used for PICO. We also introduce our new design procedure. We compare results from these three design procedures in Section 3.

The CC design approach is "Orthogonal Array Latin Hypercube Sampling," or OALHS (Tang 1993; Morris & Mitchell 1995; Leary et al. 2003; Carnell 2009). It begins with defining the boundaries of a hypercube in a scaled parameter space—all parameter dimensions rescaled to fit in an interval between 0 and 1. A grid in this hypercube is then defined. In a two-dimensional parameter space the n_d design points are assigned to the grid in such a way that each design point is the only occupant of its row and the only occupant of its column. In a three-dimensional space, the occupation of site {i, j, k} means that sites {l, j, k} are empty for all l ≠ i, sites {i, l, k} are empty for all l ≠ j, and sites {i, j, l} are empty for all l ≠ k. The generalization to n-dimensional spaces is straightforward. The final step is to then set the exact location of each design point within its grid cell. This is done using various criteria that generally aim to maximize the distance between points in any lower-dimensional projection. The fact that the points are well spread out, in many different ways of quantifying "spread out," improves interpolation to non-design points.

The design for PICO is constructed from samples from a prior probability distribution for the parameters and hence we call this design procedure PS for "Prior Sampling." In particular, the PICO design is a draw of samples from a converged MCMC run. The design is motivated by the fact that the samples are drawn from exactly the region of the parameter space where we desire accurate emulation (i.e., the region with non-negligible prior probability).

The parameter bounds used in our OALHS design are given in Table 2. The bounds have a width of eight times the marginal Fisher matrix errors for the WMAP5 noise model, centered on the fiducial point given in Column 2 of Table 1. Given these bounds and the number of design points, n_d, we draw a realization of an OALHS using the maximinLHS function of the R lhs package (Carnell 2009), which maximizes the minimum distance between the design points within the Latin hypercube. All two-dimensional projections of a 100-point OALHS design are shown in the upper triangle of panels in Figure 3.

Zoom In Zoom Out Reset image size

Both the OALHS and PS design procedures have their advantages. The OALHS design fills the parameter space in an efficient way, providing well-spaced anchor points for interpolation. The PS design, on the other hand, concentrates design points only in the region of parameter space where we are actually interested in evaluating the emulator. We now present a design procedure that captures the advantages of both.

If p_θ is the number of input parameters to our simulation, we can compute a Fisher information matrix (based on the data error distribution determining the parameter prior) from the summary statistic of the simulations we wish to emulate using 2p_θ + 1 additional simulation runs; that is, via p_θ two-point numerical derivatives of the summary statistic and one fiducial point used to calculate the summary statistic covariance model. We use this Fisher matrix to rotate the input parameter axes to a less-correlated set of input parameters. We then build an OALHS design in this de-correlated parameter space, but reject all samples of the OALHS design that are outside of a hypersphere of constant probability (at quadratic order) centered on the fiducial parameter location. We call this design OALHSFS for OALHS samples inside a "Fisher sphere." We also use this decorrelated space for performing the interpolation of the simulation design runs.

The two-dimensional projections of a 100-point OALHSFS design constructed from a Fisher matrix using the WMAP5 noise model are shown as the black triangles in the lower panels of Figure 3. The design points are located within the hypersphere of radius 4 in the decorrelated parameter space (note that each parameter has unit variance in this space). We chose a radius of 4σ because this encloses over 98% of the volume of a six-dimensional Gaussian distribution, which we judged to be sufficient for our emulator. (Note that this is covering 98% of the volume of the prior distribution. The fraction of the posterior distribution covered will be much larger.) The points of a 100-point PS design selected as random samples from an MCMC chain using the CMB power spectrum with the same WMAP5 noise model are overplotted as red crosses for comparison. It is readily apparent that the OALHSFS design points are confined to the region of high prior probability as defined by the MCMC samples.

By choosing design points inside a hypersphere in a de-correlated parameter space, we can achieve large reductions in the volume of the space that must be covered by the design points as the dimensionality of the parameter space increases. In three dimensions, the volume of a sphere embedded in a cube, such that the sphere's diameter is equal to the length of a side of the cube, is nearly half the volume of the cube. In six dimensions, the hypersphere has a volume nearly 12 times smaller than that of the surrounding hypercube. In principle then, if the interpolation error in our emulator is largely determined by the mean distance between design points, we should be able to achieve similar interpolation precision using 12 times fewer design points in an OALHSFS design as in an OALHS design in a six-dimensional parameter space. Or, for a fixed number of design points, the OALHSFS design should have smaller interpolation errors than the OALHS design.

As shown in Figure 3, the PS design gives similar volume savings over the OALHS design as does the OALHSFS. However, the PS designs do not consistently cover the design space and can lead to very large interpolation errors in regions where the design points are sparse.

The Fisher matrix approximates the shape of the likelihood about its peak at quadratic order in the cosmological parameters. Our OALHSFS design will fail to encompass a region of constant probability if higher order corrections to the shape of the likelihood are significant. The PS design does not suffer from this problem. In practice, it will be possible to measure the accuracy of the OALHSFS probability contours by comparing with a PS design as we have done in the lower panels of Figure 3. If necessary, the OALHSFS design can be adjusted to allow for different ranges along the different axes in the "de-correlated" parameter space. That is some axes might extend for, e.g., ±4σ while some might extend for ±12σ, etc. (where σ is the marginal error determined from the Fisher matrix).

In addition to the prior parameter ranges, the parameter ranges along each axis of the simulation design could also be informed by the curvature of the mode amplitude surfaces to be interpolated. In the case that a design axis has a range much shorter than the response surface curvature scale along a given coordinate axis the design will provide an estimate of the curvature scale. This could result in, e.g., the response surface being modeled as less smooth than it really is. Information about the curvature of the response surfaces is contained in the second derivatives of the power spectra with respect to the cosmological parameters (while the Fisher matrix uses the first derivatives of the power spectra). The range of the design axes could then be increased in directions where the curvature is very small (e.g., using the mode amplitude with the smallest curvature to set the axis length). On the other hand, if the response surface is known to have an exceptionally high curvature along a given axis, the density of design points could be increased to accurately sample the high curvature. Because the emulators for our case study and Emu CMB work to an acceptable precision we have not explored these ideas further.

It is always possible that when inferring parameter constraints from a data set via MCMC, the MCMC chain will move some steps out of the region of parameter space covered by the emulator design. If the MCMC chain spends very little time in the region outside the design, it may possible to fall back on evaluating the actual Boltzmann code for analyzing the CMB. For applications with more expensive simulations (such as N-body simulations), the emulator's predictions can be extrapolated to get an initial MCMC chain that can then be used to inform the construction of a new design.

To reduce the dimensionality of the function to be interpolated, both the CC and PICO emulators perform a principle component decomposition of the n_y × n_d matrix of simulation design runs yⁱ_ℓ ≡ ℓ(ℓ + 1)/(2π) Cⁱ_ℓ where i = 1, n_d, and ℓ takes n_y values in the range (2, ℓ_max). The principal components (PCs) that contribute least to the variance over the design are discarded. We have found that keeping the 20 most important modes is more than sufficient for our purposes, i.e., the errors made from neglecting the other min(n_y, n_d) modes are insignificant. We quantify the errors from the truncated mode decomposition in two ways. First, we have verified that using 20 PC modes allows us to reconstruct the power spectra at all design points in several OALHS designs to within 0.1%. Second, we used the Fisher matrix to estimate the cosmological parameter biases that would be induced from the systematic errors in the reconstructed power spectra (see Huterer & Takada 2005 for details on the Fisher matrix methods). For p_μ = 20, the distribution of parameter biases had first and third quartiles <10% for all six of our cosmological parameters. This means that in the case of perfect interpolation of the power spectra, the inferred parameter constraints using our emulator would have biases less than 10% of the 1σ marginal errors for both the OALHS and OALHSFS designs.

Following Heitmann et al. (2006), Habib et al. (2007), and Schneider et al. (2008), we first subtract the (ℓ-dependent) mean of the design runs from each design power spectrum and then scale the result by a single number so that the combined entries of our n_y × n_d matrix of design runs have variance one. Denoting this centered and scaled matrix η, we then perform a singular value decomposition, η = UBV^T, where U has dimension n_y × p (p ≡ min(n_y, n_d)) with $\mathbf {U}^{T}\mathbf {U}=\mathbb {I}_{p}$, V has dimension n_d × p with $\mathbf {V}^{T}\mathbf {V}=\mathbb {I}_{p}$, $\mathbf {V}\mathbf {V}^{T}=\mathbb {I}_{n_{d}}$, and B (p × p) is a diagonal matrix of singular values. Finally we define the matrix of basis vectors, ${\bPhi }\equiv \frac{1}{\sqrt{n_{d}}}\mathbf {U}\mathbf {B}$ and $\boldmath {w}\equiv \sqrt{n_{d}}\mathbf {V}^{T}$ so that $\frac{1}{n_{d}}\boldmath {w}^{T}\boldmath {w}=\mathbb {I}_{n_{d}}$.

Retaining only the p_μ ⩽ p most significant columns of ${\bPhi }$ we can write the scaled power spectrum as a sum over p_μ modes,

$$\begin{equation} \eta _{\ell }({\theta }) = \sum _{\mu =1}^{p_{\mu }} w_{\mu }({\theta }) \Phi _{\ell \mu }, \end{equation} \tag{ 4 }$$

where Φ_ℓμ is the entry of the matrix ${\bPhi }$ in row ℓ and column μ and $w_{\mu }(\boldmath {\theta })$ is the μth (parameter-dependent) mode amplitude.

In Figure 4, we show the first four PC basis functions for the power spectra in the three different design types, PS, OALHS, and OALHSFS. Because of the very different volumes of parameter space covered, application of the PC algorithm to the design runs yields very different basis functions for the OALHS design versus the PS and OALHSFS designs. This is a crucial difference because the basis functions must accurately describe the design spectra at all multipoles in order to meet the accuracy targets over the large dynamic range we are simulating. Note that modes 1 and 2 derived from the PS and OALHSFS designs have several oscillations, but are much smoother when derived from the OALHS design. This difference can be seen later in the interpolation errors from the different designs (Figure 7) where the OALHS design has larger, and more oscillatory, errors at large multipoles indicating that the missing oscillatory structure in the OALHS design basis functions leads to missing structure in the reconstructed power spectra. In Appendix A, we address the choice of basis functions again when we extend the emulator for the CMB temperature power spectrum to ℓ = 5000.

Zoom In Zoom Out Reset image size

We compare two methods for interpolating the mode amplitudes $w_{\mu }(\boldmath {\theta })$ between design points in the cosmological parameter space.

For the first interpolation method we again follow the CC framework and model $w_{\mu }(\boldmath {\theta })$ as independent GPs for each value of μ. We refer the reader to Habib et al. (2007) and Heitmann et al. (2009) for details about the construction of the GP model for the mode amplitudes. In summary, for a fixed set of points $\boldmath {\theta }_i$ (i = 1, ..., n_d), the amplitudes w_μ(θ_i) are modeled to have a multivariate Gaussian distribution with mean zero and a parameterized covariance structure that determines the smoothness of the interpolation. The covariance parameters are calibrated from a likelihood model conditioned on the simulation design runs. Because the GP framework infers distributions for the covariance parameters rather than best fit values, the error in the GP fit and interpolation can be propagated through to the final cosmological parameter constraints.

Our second interpolation method is to fit a p_θ-dimensional linear or quadratic polynomial to the design runs. Compared to the GP model, this has the advantage of being very fast to compute, but the disadvantage of not propagating any interpolation error. For small numbers of design points, the order of the polynomial is limited (in order to obtain a well-conditioned fit) and this can further limit the accuracy of the polynomial interpolation. We therefore expect GP interpolation to outperform global polynomial fits when the curvature of the mode amplitude surfaces (in any parameter direction) exceeds the maximum polynomial order that can be fit with n_d design points. A caveat is that GP models with our chosen covariance will also yield unacceptably large interpolation errors for extremely rapidly varying surfaces unless the number of design points is significantly increased. So there is a general limit on the smoothness of the mode amplitudes as functions of cosmological parameters for such high-dimensional parameter spaces and limits on the practical number of design runs.

Figure 5 shows the first and fourth mode amplitudes in the PC decomposition of the power spectra in the OALHS 100-point design as a function of n_s with all other cosmological parameters fixed at their central values. The range of n_s covers the marginal 4σ error as derived from the Fisher matrix with the WMAP5 noise model. The first mode amplitude in this case is a very smooth function of n_s, while the less significant modes have a higher degree of curvature. (Note that the size of the residuals in the mode amplitudes does not obviously translate into power spectrum errors because of the centering and scaling of the power spectra done before performing the mode decomposition.)

Zoom In Zoom Out Reset image size

For our case study we have found the step of calibrating the GP parameters to have intensive computational and human-labor requirements. With the GP covariance parameterization used in the CC framework there are p_μ variance parameters and p_μ × p_θ correlation parameters. To calibrate the GP model we must then find the maximum of the likelihood in this high-dimensional parameter space. We used MCMC to find the vicinity of the GP likelihood peak, but found the most difficult step to be initializing the chains in regions of non-negligible probability. A simulated annealing algorithm with a constant cooling schedule turned out to be quite effective in getting our GP calibration chains started. Because the calibration of the GP model can be a large computational obstacle in constructing an emulator we highly recommend the use of a simulated annealing or similar algorithm. We do not, however, know of any algorithm that can replace the human labor of assessing the convergence of the multiple MCMC chains required to optimize the GP model. We demonstrate the performance of our simulated annealing algorithm for calibrating the GP interpolation parameters in Figure 6.

Zoom In Zoom Out Reset image size

We also found the choice of the proposal distribution for the GP correlation parameters in the MCMC algorithm to be crucial to efficient calibration of the GP model. We used independent uniform distributions for each correlation parameter, but initially allowed the width of the uniform proposal to be adjusted according to the variance in the previous 1000 steps in the chain. Because this breaks the detailed balance requirement for convergence of the MCMC chain, we fix the widths of the proposal distribution after a preset number of chain steps (of order 10⁵).

Because the evaluation of the GP likelihood at each chain step requires inversion of an n_d × n_d matrix, the GP interpolation method becomes unpractical as n_d becomes very large. This is unfortunate because the interpolation accuracy should improve as n_d increases! We therefore expect GP models to provide a preferable interpolation method for small n_d and polynomial interpolation to be preferred for large n_d. For our case study, we actually find that polynomial interpolation provides smaller interpolation errors for a wide range of n_d. We demonstrate the performance of polynomial interpolation for a different case study matching that of Heitmann et al. (2009) in Appendix B.

Our C++ code for calibrating GP models for interpolation (built with the Scythe Statistical library; Pemstein et al. 2007) is available for download at http://www.emucmb.info.

In Figure 7, we show the fractional interpolation error for 100 test spectra using the three design methods under consideration. The test spectra were computed at points in our six-dimensional cosmological parameter space drawn from the posterior distribution for our WMAP5 noise model. Each interpolation method was trained on 100 design points. (Note that we always show the expected value of the emulator; we do not compute the random draws that are part of a complete GP prediction.) At these test points, the OALHSFS design has interpolation errors approximately two times smaller than the OALHS design at all multipoles and for both interpolation methods. Using quadratic polynomial interpolation, the PS design has comparable interpolation errors to the OALHSFS design, but has very poor performance using GP interpolation. As we mentioned in Section 2.2, the OALHSFS design is built in a parameter-space volume that is ∼1/12 the size of the OALHS volume. If the increased density of simulation points is the dominant difference between the OALHSFS and OALHS designs, then Figure 7 gives some indication of how the interpolation accuracy improves with the increased density. However the interpretation is not straightforward as the shape of the designs (spherical versus cubic) may also be important. That is, the design points outside the Fisher sphere in the OALHS design could, in principle, be helpful in reducing the interpolation accuracy, but we have not pursued the separation of these effects further. The main result that the OALHSFS design gives improved performance is sufficient for our purposes. We show similar results, but for n_d = 50 in Figure 8. The interpolation errors using quadratic polynomials have similar relative performance between the three design types as in Figure 7. The GP interpolation results, however, now show smaller errors using the PS design than using the OALHS design. The OALHSFS design using GP interpolation generally gives smaller errors than the OALHS design, but there are a few test points where the OALHSFS design gives worse performance for several multipoles. These results are consistent with our expectation that at a fixed number of design points n_d the reduction in the mean distance between design points in the OALHSFS and PS designs (see Figure 3) should improve the interpolation accuracy. Indeed we do see an improvement by a factor of ∼2 in the interpolation accuracy. However, the poor performance of the PS design with 100 design points and using GP interpolation demonstrates that the PS design is a less reliable algorithm than OALHS-based designs for building an emulator using GP interpolation. This could be because the PS design leaves regions of parameter space very sparsely covered compared to the OALHS designs (especially in higher dimensions).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Next, we investigate how the accuracy gains of OALHSFS over OALHS at fixed n_d translate into a reduction in n_d at fixed interpolation accuracy. In Figure 9, we show the quartiles and maxima of the average fractional interpolation error in a band from ℓ = 1500 to ℓ = 2000 (that is, we averaged the absolute value of the error over ℓ = 1500 to 2000) versus n_d for 100 test points taken from a PS design. We plot the interpolation errors for both the OALHS and OALHSFS designs using GP interpolation. Comparing the solid and dashed lines at fixed interpolation accuracy (i.e., along a horizontal axis) we see the OALHSFS design requires ∼25%–70% fewer design points than the OALHS design (with larger savings at low n_d). One possible explanation for this trend is that some of the design points in the OALHS design are not significantly contributing to the reconstruction of the mode amplitude surfaces in the (high-probability) region where the test points are actually located. When n_d is small, every design point should be important for reducing the interpolation errors, so wasting even one point in the OALHS design would degrade the emulator. This demonstrates a main point of the paper that the design space volume reduction offered by the OALHSFS design is significant in minimizing n_d to achieve a fixed error tolerance.

Zoom In Zoom Out Reset image size

In order to understand when the errors in the truncated mode decomposition of the power spectra are significant relevant to the interpolation errors, we plot the fractional errors in the reconstructed power spectra as functions of p_μ at each design point (i.e., without any interpolation) in Figure 10. The reconstruction errors are rapidly decreasing functions of the number of retained PC modes, as expected. For the choice used throughout this paper of p_μ = 20 for n_d ⩾ 50, the errors from the truncated mode decomposition are adequate for our error tolerance of 10⁻³ for both the OALHS and OALHSFS designs. However, from Figure 9 we can see that with p_μ = 20 the mode reconstruction errors are non-negligible for the n_d = 100, 200 designs, which would explain why the total fractional power spectrum errors (mode reconstruction + interpolation) do not steadily decrease for n_d > 100.

Zoom In Zoom Out Reset image size

In Figure 11, we compare the mean interpolation errors in the band 1500 < ℓ < 2000 using the GP model and a quadratic polynomial fit to the OALHSFS and PS designs as functions of n_d. In all of the designs considered polynomial interpolation provides smaller interpolation errors than GP interpolation for small n_d. For the OALHSFS design, the interpolation methods give comparable errors for large n_d. The PS designs have sporadic error ranges, indicating the large variation in the covering of the design space test points with different PS design realizations.

Zoom In Zoom Out Reset image size

The (first- or second-order) polynomial interpolation performs so well because the response surfaces of each mode amplitude are so smooth as functions of the six active cosmological parameters. The GP models also require smooth response surfaces, but we expect the low-order polynomial fits to be worse than the GP fits for a given correlation length in the response surface as the correlation length is decreased. This is because the covariance structure in the GP model has more flexibility than a (global) second-order polynomial. Therefore, for some applications the fast-to-compute polynomial interpolation will no longer suffice. The relative performance of different interpolation methods will always be problem specific. However, we explore designs for the matter power spectrum in Appendix B and show that the polynomial interpolation is again comparable to the GP interpolation demonstrated in Heitmann et al. (2009) for similar designs for most of the dynamic range in the simulated power spectra (except for the smallest scales). So, for both the CMB and matter power spectrum, an emulator constructed using an OALHSFS design and quadratic polynomial interpolation has similar performance to an emulator using an OALHS design and GP interpolation. As the number of design points increases, both interpolation methods should converge to arbitrary accuracy as higher order polynomial fits become well conditioned.

We use an OALHSFS design with n_d = 100 covering the WMAP5 4σ error bounds as determined by our Fisher matrix described in Section 2.1. The design runs for Emu CMB use the 2009 October release of CAMB (different from the case study in the main body of the paper) with the "curved correlation function" lensing implementation including nonlinear corrections calculated by HALOFIT (Challinor & Lewis 2005). We again use accuracy settings accuracy_boost = 4 and l_accuracy_boost = 4, but build the emulator from the power spectra as output by CAMB at every multipole rather than extracting only the uninterpolated power spectra values. We set l_sample_boost = 4 to determine the number of multipole values where the power spectra are actually computed. See http://www.emucmb.info for the complete CAMB parameter files.

Also in contrast to the main body of the paper, we perform the dimension reduction of the design runs on yⁱ_ℓ ≡ log(ℓ(ℓ + 1)/(2π) Cⁱ_ℓ) (i = 1, n_d). Performing a PC decomposition on the logarithm of the power spectra yields basis functions that better describe the lensing variations at ℓ ≳ 3000 across the design points as shown in Figure 12.

Zoom In Zoom Out Reset image size

As in the main text, we fit global second-order polynomials to each of the first 20 mode amplitudes to perform the interpolation over the design space. Using only the first 20 PCs allows reconstruction of the power spectra at the design points to less than 0.1% accuracy for ℓ ≲ 100 and to less than 0.01% accuracy for ℓ ≳ 100.

We validated the emulator by running CAMB at an additional 100 points in a new realization of an OALHSFS design. (Note the validation plots in the main body of the paper use points from a PS design.) The interpolation errors for the spectra at these 100 test points are shown in Figure 13. For all test points, the emulator is accurate to less than 0.6% out to ℓ = 5238. The errors are less than 0.2% for ℓ < 3000, which is consistent with Figure 7.

Zoom In Zoom Out Reset image size

Note that while Emu CMB uses significantly fewer training runs of CAMB than PICO, unlike PICO this Emu CMB example does not predict polarization power spectra. We have verified with simple tests that interpolation of the low-ℓ parts of the EE and TE CMB power spectra is more challenging than the TT spectra demonstrated here. In addition, many of the training runs for PICO were required to include cosmological models with non-zero spatial curvature, which we have not considered.¹²