Optimising Boltzmann codes for the Planck era
Abstract
High precision measurements of the Cosmic Microwave Background (CMB) anisotropies, as can be expected from the planck satellite, will require highaccuracy theoretical predictions as well. One possible source of theoretical uncertainty is the numerical error in the output of the Boltzmann codes used to calculate angular power spectra. In this work, we carry out an extensive study of the numerical accuracy of the public Boltzmann code CAMB, and identify a set of parameters which determine the error of its output. We show that at the current default settings, the cosmological parameters extracted from data of future experiments like Planck can be biased by several tenths of a standard deviation for the six parameters of the standard CDM model, and potentially more seriously for extended models. We perform an optimisation procedure that leads the code to achieve sufficient precision while at the same time keeping the computation time within reasonable limits. Our conclusion is that the contribution of numerical errors to the theoretical uncertainty of model predictions is well under control – the main challenges for more accurate calculations of CMB spectra will be of an astrophysical nature instead.
LAPTH1317/09
1 Introduction
Any meaningful quantitative analysis of experimental data is based on a comparison with the predictions of a theoretical model, and the analysis of Cosmic Microwave Background (CMB) anisotropies is no exception to this rule. The amount of information to be gained from observations is limited not only by experimental uncertainties (such as detector noise, foregrounds, etc.), but also by the ability to accurately predict observable quantities from a given theory. There are various sources of theoretical uncertainties. Some, such as cosmic variance, are endemic to the problem, and unavoidable. Others are based on insufficient theoretical understanding of the complex processes involved [1] (examples include the physics of recombination [2], reionisation [3], and contributions due to the SunyaevZel’dovich effects). Additionally, the analysis may be compromised by inaccuracies in the numerical programmes used to calculate the anisotropy angular power spectra [4].
With increasingly sophisticated experiments, and the contribution from experimental errors becoming less and less important, the relative significance of theoretical errors increases. In fact, for the planck satellite [5], the signal uncertainty of the temperature anisotropies will be dominated by cosmic variance instead of noise over a wide range of scales up to multipoles of .
Ignoring any type of uncertainty can lead to biased estimates of parameter values, and, in the worst case, a wrong physical interpretation of the data. It is therefore imperative to either make sure that the errors are small enough to be negligible, or, if that should not be the case, to devise an appropriate strategy to deal with the problem.
In the present work, we shall consider the numerical accuracy of the Boltzmann codes employed to calculate the angular anisotropy spectra for given input values of cosmological parameters. The first public Boltzmann code was released more than a decade ago [6], and to date, there are several other such programmes freely available for download [7, 8, 9]. The output of all these codes is necessarily inaccurate to some extent, due to the use of semianalytical approximations as well as artifacts of the numerical implementation, such as finite integration steps or the need to interpolate. These effects can be parameterised by a set of accuracy parameters, whose settings determine the accuracy of the output, but also the computation time. Here, we will focus our analysis on the CAMB code by Lewis, Challinor and Lasenby [8] in order to avoid possible systematic effects caused by differences between codes – a comparison of different (more or less) independent Boltzmann codes was performed by Seljak et al. [4], who found an excellent qualitative agreement.^{1}^{1}1We verified that after updating various parts of CMBfast (values of physical constants, recombination code, etc.), the output of CAMB and CMBfast agree sufficiently well. We extend their line of reasoning and present a detailed analysis of the potential effects of numerical inaccuracies on parameter estimates from planck data. Additionally, we optimise the accuracy settings of CAMB to find an ideal balance between precision and execution time.
We proceed in Section 2 by defining an appropriate measure of accuracy, identifying the relevant parameters which affect the accuracy of the output spectra and constructing a set of CMB reference spectra. In Section 3, we will describe our optimisation procedure and present a recommended set of accuracy parameters for CAMB, followed by an analysis of the expected potential bias on the cosmological parameters of the vanilla model caused by numerical inaccuracies in Section 4. The impatient reader may prefer to skip directly to our conclusions in Section 5.
2 Reference spectra
2.1 The fiducial model
In principle, the impact of individual accuracy parameters and the number of relevant parameters will depend on the underlying cosmological model assumed, and on the values of the cosmological parameters. In this analysis, we will stick to the 6parameter CDM“vanilla”model and we limit ourselves to a point in the space of cosmological parameters that lies close to the best fit to the WMAP 5year data [10], see Table 1. We neglect the effects of weak gravitational lensing on the CMB spectra [11] (which will of course have to be taken into account in an analysis of real planck data), and ignore signatures of nonminimal models, like massive neutrinos, tensor modes, spatial curvature, etc., and defer their treatment to future work.
Parameter  Fiducial Value  Prior Range  

Dark matter density  0.10976  
Baryon density  0.02303  
Hubble parameter  0.7  
Optical depth to reionisation  0.09  
Scalar spectral index  0.96  
Amplitude of scalar spectrum @ Mpc  3.135 
2.2 Measuring accuracy
In order to quantify the accuracy of the Boltzmann code output , we require two things: a reference point to compare with, and a measure of accuracy.
The reference spectra would ideally be the exact prediction of the theory. In practice however, we have to make do with getting close enough to these ideal spectra. We will return to this issue and describe the construction of in 2.4.
To measure accuracy, one might be tempted at first glance to look at the relative difference of the spectra for each . However, this approach does not properly take into account the fact that the accuracy requirements are dependent on (due to cosmic variance and experimental errors). Additionally, the spectra are merely an intermediate step in the inference process. In the end, we are interested in possible biases on cosmological parameters, not the accuracy of the spectra. A more meaningful measure of deviation from the reference spectra is the effective , which can be obtained by using the reference spectra to generate a fiducial data set, taking into account the experimental errors of an experiment (or the projected errors in case of future experiments), and “fitting” to these data.
More formally, the effective is related to the likelihood and defined by
(2.1) 
Here, is the theoretical covariance matrix, and its entries are taken to be the sum of the signal and noise power spectra:
(2.2) 
where the index runs over temperature () and polarisation (), and for a fiducial data set, we can take the data covariance matrix to be equal to , i.e., the theoretical covariance matrix evaluated for the fiducial values of the cosmological parameters.
We take the noise to be isotropic and Gaussian; the noise power spectrum is related to the experimental parameters of Table 2 through
(2.3) 
For more details see Refs. [12, 13]. It should be noted that the normalisation of Eq. (2.1) is chosen such that the total is zero when the output spectra exactly match the reference spectra used to construct the fiducial data.
We generate the fiducial data set of , , and spectra up to using the code of Perotto et al. [12]. For simplicity we ignore the effects of incomplete sky coverage due to masking the galaxy and point sources, as well as anisotropic noise [14]. To evaluate the accuracy of CAMB in view of the expected data from planck, we assume 14 months of integrated observations in the 70 GHz channel of LFI and the 100 and 143 GHz channels of the HFI instrument; their specifications are taken from the planck blue book [5] and listed in Table 2.
2.2.1 Interpretation of the measure
As can be seen in Eq. 2.1, is directly related to the likelihood, which, along with a choice of prior probability densities on all cosmological parameters, leads to the posterior probability density, from which constraints on parameters are eventually derived. Assuming flat priors on the parameters and a multivariate Gaussian likelihood function, for a given numerical error , the bias on any cosmological parameter cannot exceed standard deviations in the worst case (i.e., when the error in the angular power spectra can be exactly offset by changing one of the cosmological parameters). On the other hand, if the error had no degeneracy with any cosmological parameter, the effect of a nonzero would be just a constant offset to the likelihood function, which would not have any effect on parameter inference. In realistic cases, the expected bias would lie somewhere in between. For the parameters of the vanilla model, we will provide an estimate of the bias in Section 4.
/GHz  /K  /K  

70  14.0’  12.8  18.3 
100  9.5’  6.8  10.9 
143  7.1’  6.0  11.4 
2.3 The accuracy parameters
The numerical accuracy of a Boltzmann code’s output depends on many factors, ranging from the use of analytical approximations to the sampling various intermediate quantities in the calculation. These sources of numerical error can be quantified in terms of accuracy parameters, e.g., the number of samples used for interpolating a particular quantity. An increase in accuracy will generally be accompanied by a longer computation time and possibly higher requirements on the available computer memory.
We use the June 2008 version of CAMB^{2}^{2}2http://www.camb.info as a starting point of our analysis. The unmodified version of CAMB comes with a set of three continuously adjustable accuracy parameters:^{3}^{3}3There are also a few logical switches that are relevant here; in our analysis we kept them fixed to accurate_polarization=T, accurate_reionization=T and do_late_rad_truncation=T.

l_sample_boost: determines for which values of the are actually calculated (the rest are interpolated).

l_accuracy_boost: determines the multipole at which the Boltzmann hierarchy for photons, neutrinos, etc., is cut off.

accuracy_boost: affects the setting of various time steps, samplings, etc.
The latter two parameters affect several settings at once, so in the interest of optimising the performance of CAMB, we split them up into their constituents and treat them separately. Apart from the settings governed by these three parameters, we identified a few other quantities which can affect the accuracy of the results, and should be taken into account when optimising the code. Altogether, we consider a set of 19 accuracy parameters in our analysis, listed in Table 3.
All parameters are defined in such a way that setting them a value of 1 reproduces the results of the unmodified version of CAMB, and larger values correspond to better accuracy. The constituents of the old l_accuracy_boost and accuracy_boost parameters are taken to multiply the old l_accuracy_boost and accuracy_boost parameters (e.g., setting produces the same effect as setting ). We modified the routine that determines for which values of the are calculated: our parameter new_l_sample_boost is defined to be the square root of the old l_sample_boost; for , all are calculated and there will be no interpolation of the final spectrum.
Parameter  Corresponding old accuracy parameter  Description and comments 
new_l_sample_boost  l_sample_boost  sampling of s 
lmaxg  l_accuracy_boost  Boltzmann hierarchy cutoff for photons 
lmaxnr  Boltzmann hierarchy cutoff for massless neutrinos  
int_tol  tolerance parameter for integration routines  
ri_timestep  time step during reionisation  
rec_timestep  time step during recombination  
rec_timestep2  time step between recombination and reionisation  
rad_trunc  truncation of photon hierarchy during matter domination  
dec_start  accuracy_boost  starting time of decoupling 
int_sample  samples for integration over source function  
source_dk  sampling of source function  
source_kmin  minimum value of to calculate source function for  
int_xlmax1  starting time for source function integration  
tc_largek  switch off tight coupling later for large  
tc_ep0  –  tight coupling switch 
bess_sampling  –  sampling of spherical Bessel functions 
bess_xlimmin  –  approximate for small , if 
bess_xlimfrac  –  approximate for large , if 
ketamax  –  maximum value of 
2.4 Constructing a reference data set
To quantify absolute accuracy we require a reference data set, as discussed above. Its construction is naturally tied to choosing the accuracy parameters in such a way that increasing them further would not have any appreciable effect. However, by arbitrarily increasing all parameters to “large” values at the same time, one would run into the limits of the hardware, particularly the available memory. For the purpose of finding suitable values for generating the reference spectra, we therefore analyse the parameters one by one, keeping all other parameters fixed. For each parameter, we generate a fiducial reference data set with that parameter set to a high value, all other parameters kept at a value of 2. Varying this parameter and calculating the reveals its impact on overall accuracy, allows us to find a suitable setting for the reference spectra and lets us estimate the remaining error.
We show the results of the single parameter scans in Figure 12. From Figure 1 we can see that certain parameters (e.g., lmaxg, lmaxnr) are very wellbehaved and reach a plateau well below the fiducial value. For these parameters it is reasonable to assume that picking even higher values would not have any appreciable effect on accuracy. A number of other parameters do not fully converge, and exhibit a steplike behaviour before reaching the fiducial value (e.g., int_sample, source_dk). It is likely that increasing their value beyond our fiducial maximum value would have an effect on the spectra. However, the graphs in Figure 1 nonetheless provide an order of magnitude estimate of the remaining error and one can use them as a guide to finding reasonable settings for the construction of the reference spectrum. Finally, the parameter source_kmin does not seem to converge at all and exhibits unstable behaviour, though its effect on overall accuracy is negligible.
The parameter settings we chose for the reference spectrum are given in Table 4. The dominant contribution to any residual error of the reference spectrum will come from the parameter displaying the worst convergence – ketamax. Unfortunately, this parameter also has a strong impact on the computation time (see Figure 2), and memory requirements, precluding us from choosing a higher setting. We estimate the reference spectrum to be accurate to of order .
Parameter  Reference setting  Setting 1  Setting 2  Setting 3 

new_l_sample_boost  6  1.77  1.60  2.00 
lmaxg  6  7.39  2.02  2.05 
lmaxnr  6  1.78  2.11  5.36 
int_tol  6  2.07  1.53  5.86 
ri_timestep  2  0.99  0.71  2.87 
rec_timestep  3  0.50  0.50  0.75 
rec_timestep2  10  1.85  1.14  2.44 
rad_trunc  4  1.81  1.76  2.58 
dec_start  100  2.35  35.76  5.05 
int_sample  6  4.21  4.88  3.65 
source_dk  4  3.10  2.84  2.72 
source_kmin  1  5.86  2.50  5.23 
int_xlmax1  4  1.20  1.00  2.16 
tc_largek  5  2.44  2.33  1.90 
tc_ep0  5  4.25  3.30  6.32 
bess_sampling  3  3.02  2.53  4.00 
bess_xlimmin  2  1.14  0.74  3.61 
bess_xlimfrac  2  4.19  0.90  1.07 
ketamax  3  1.32  0.96  0.56 
3 Optimising performance
Having constructed the reference spectra, we can now proceed to searching settings of the accuracy parameters which yield a result that lies as close as possible to the reference spectra, within a reasonable time of computation. To this end, we use a modified version of the Markov chain Monte Carlo code CosmoMC [15]. Major modifications include:

we use the fiducial data set constructed from the reference spectra;

we vary the accuracy parameters instead of the cosmological parameters (which are kept fixed at their fiducial values), taking top hat priors on all accuracy parameters (with lower limits at a value of 0.5 and upper limits large enough to not influence the results);

instead of sampling from the usual posterior (which is proportional to ) itself, we sample a function , defined by
(3.1) This function was chosen such that areas of parameter space leading to too long computation times are avoided, and that areas of parameter space giving are better sampled.
We generated sample settings in this way; a scatter plot of the samples in the plane is presented in Figure 3, illustrating the strong correlation between accuracy and computation time. Which settings to use is a somewhat subjective decision and should be taken with one’s available computing power in mind. We list two sample settings in Table 4, taken from the samples of our Monte Carlo run: “Setting 1” corresponds to the best accuracy under the constraint that the computing time be less than 60 seconds, “Setting 2” gives the best accuracy for , and “Setting 3” the best accuracy for (which is the computation time for the unmodified version of CAMB run at accuracy_level = 2). The performance of these settings is contrasted to CAMB’s default (accuracy_level = 1), and a highaccuracy setting of an unmodified version of CAMB (accuracy_level = 2) in Table 5.
Note that at the default settings, with a difference of to the reference spectra, parameter estimates could be biased by more than two standard deviations, in the worst case. For both Settings 1 and 2, on the other hand, the maximum possible bias would be less than 0.1 standard deviations ( standard deviations for Setting 3), assuming Gaussian posterior distributions.
The results of the MCMC search confirm the tendencies observed in the single parameter scans of Figs. 1 and 2 regarding the impact of individual parameters on accuracy and speed; we find no evidence for significant crosscorrelations between accuracy parameters.
Setting  

1  59  
2  27  
3  17.4  
5.8  2.7  
0.5  17.6 
4 Estimating the bias
As mentioned above, the for a given accuracy settings can be used to estimate the worst case bias on cosmological parameters. Nevertheless, we are also interested in how large a bias effect we can expect in the example case of the vanilla model. To this end, we performed an actual parameter estimation exercise using CosmoMC, for three cases:

Using the same CAMB accuracy settings for generating the fiducial data as for the inference process – mimicking an analysis free of numerical errors.

Using the reference data set and running CAMB at accuracy settings “2” from Table 4.

Using the reference data set and running CAMB at its default settings.
We generate 16 Markov chains, making sure the GelmanRubin convergence parameter [16] is smaller than for all parameters considered. The results are plotted in Figure 4: as expected, there is no discernible bias between the results obtained with optimised settings and an errorfree analysis. The expected bias for the default settings is rather mild, not exceeding a few tenths of the posterior standard deviations for the respective parameters, with (40%) and (43%) being the most affected.
5 Conclusions
To meet the challenge posed by ultraprecise future CMB experiments, the output of the Boltzmann codes used to calculate theoretical predictions for the CMB anisotropy spectra will have to meet stringent requirements in terms of numerical accuracy – a goal that comes at the cost of accordingly higher demand in computational resources. This tendency is partially mitigated by the increase in available computing power. Additionally, the use of efficient interpolation algorithms [17, 18, 19, 20] can significantly accelerate the process of parameter estimation from future data sets. Nonetheless, such methods still rely on the input of a full Boltzmann code, whose inherent numerical errors will be propagated to the interpolation codes.
We have performed a detailed analysis of the numerical accuracy of CAMB, provided an estimate of the residual numerical error of the output, and evaluated the possible impact on parameter estimates from planck data. For the present default accuracy settings of CAMB and simulated planck data, numerical errors can lead to a slight bias on estimates of the six free parameters of the CDM model; though we cannot exclude the possibility of a more serious bias on parameters of extended models. However, by tweaking the settings of various internal parameters of CAMB it can be made sure that the bias on any parameter will not exceed 0.1 standard deviations. An example for suggested settings of the accuracy parameters is given in Table 4.
The results of this paper lead to an efficient and more accurate calculation of CMB angular power spectra, and should bring CAMB to a standard that will allow us to make the most out of upcoming planck data. We conclude that the contribution of numerical errors to the theoretical uncertainty of model predictions is well under control – the main challenges for more accurate calculations of CMB spectra will be of an astrophysical nature instead.
Acknowledgments
We thank Martin Bucher, Anthony Challinor, Loris Colombo, Fabio Finelli and Antony Lewis for discussions and helpful suggestions.
References
References
 [1] W. Hu, D. Scott, N. Sugiyama and M. J. White, Phys. Rev. D 52 (1995) 5498 [arXiv:astroph/9505043].
 [2] W. A. Fendt, J. Chluba, J. A. RubinoMartin and B. D. Wandelt, arXiv:0807.2577 [astroph].
 [3] A. Lewis, J. Weller and R. Battye, Mon. Not. Roy. Astron. Soc. 373 (2006) 561 [arXiv:astroph/0606552].
 [4] U. Seljak, N. Sugiyama, M. J. White and M. Zaldarriaga, Phys. Rev. D 68 (2003) 083507 [arXiv:astroph/0306052].
 [5] [Planck Collaboration], arXiv:astroph/0604069.
 [6] E. Bertschinger, arXiv:astroph/9506070.
 [7] U. Seljak and M. Zaldarriaga, Astrophys. J. 469 (1996) 437 [arXiv:astroph/9603033].
 [8] A. Lewis, A. Challinor and A. Lasenby, Astrophys. J. 538 (2000) 473 [arXiv:astroph/9911177].
 [9] M. Doran, JCAP 0510 (2005) 011 [arXiv:astroph/0302138].
 [10] J. Dunkley et al. [WMAP Collaboration], arXiv:0803.0586 [astroph].
 [11] A. Lewis and A. Challinor, Phys. Rept. 429 (2006) 1 [arXiv:astroph/0601594].
 [12] L. Perotto, J. Lesgourgues, S. Hannestad, H. Tu and Y. Y. Y. Wong, JCAP 0610 (2006) 013 [arXiv:astroph/0606227].
 [13] J. Hamann and Y. Y. Y. Wong, JCAP 0803 (2008) 025 [arXiv:0709.4423 [astroph]].
 [14] S. Hamimeche and A. Lewis, Phys. Rev. D 77 (2008) 103013 [arXiv:0801.0554 [astroph]].
 [15] A. Lewis and S. Bridle, Phys. Rev. D 66 (2002) 103511 [arXiv:astroph/0205436].
 [16] A. Gelman and D. B. Rubin, Statist. Sci. 7 457511
 [17] H. B. Sandvik, M. Tegmark, X. M. Wang and M. Zaldarriaga, Phys. Rev. D 69 (2004) 063005 [arXiv:astroph/0311544].
 [18] W. A. Fendt and B. D. Wandelt, Astrophys. J. 654 (2006) 2 [arXiv:astroph/0606709].
 [19] T. Auld, M. Bridges, M. P. Hobson and S. F. Gull, Mon. Not. Roy. Astron. Soc. Lett. 376 (2007) L11 [arXiv:astroph/0608174].
 [20] W. A. Fendt and B. D. Wandelt, arXiv:0712.0194 [astroph].