© ESO, 2014

1. Introduction

The primordial power spectrum encodes the physics of the early Universe and its measurement is one of the most important research areas in modern cosmology. Amongst the proposed models that describe the early Universe, inflation (Guth 1981; Linde 1982) is currently the most favoured one. In this model early perturbations are produced by quantum fluctuations during the epoch of an accelerated expansion. These perturbations then grow into the large scale structure we observe today. The simplest models of inflation predict almost purely adiabatic primordial perturbations with a near scale-invariant power spectrum. In these models the primordial power spectrum is often described in terms of a spectral index n_s and an amplitude of the perturbations A_s as $\mathit{P}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}{\mathit{A}}_{\mathrm{s}}{\left(\frac{\mathit{k}}{{\mathit{k}}_{\mathrm{p}}}\right)}^{{\mathit{n}}_{\mathrm{s}}\mathrm{-}\mathrm{1}}\mathit{,}$(1)where k_p is a pivot scale. This spectrum represents the initial conditions set at inflation. The simplest ansatz for characterising the primordial perturbations is the so-called Harrison-Zel’dovich (HZ) model, which sets n_s = 1 (Harrison 1970; Zel’dovich 1972). This is an exact scale-invariant spectrum and has been ruled out by different datasets, as will be discussed later. Instead, the near scale-invariant spectrum with n_s < 1 fits the current observations very well (e.g., Planck Collaboration XXII 2014). However, numerous models have been proposed for the generation of the perturbations, predicting deviations from the perfectly scale-invariant power spectrum. The simplest are the slow-roll inflationary models which describe the deviations through a minimal scale dependence of the spectral index of the power spectrum, the so-called running α_s, formulated as $\mathit{P}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}{\mathit{A}}_{\mathrm{s}}{\left(\frac{\mathit{k}}{{\mathit{k}}_{\mathrm{p}}}\right)}^{{\mathit{n}}_{\mathrm{s}}\mathrm{-}\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\alpha }}_{\mathrm{s}}\ln \left(\mathit{k}\mathit{/}{\mathit{k}}_{\mathrm{p}}\right)}\mathrm{·}$(2)More complex models generating deviations from scale invariance include those with features on the potential (Starobinsky 1992; Adams et al. 2001; Wang et al. 2005; Hunt & Sarkar 2004, 2007; Joy et al. 2008; Pahud et al. 2009; Lerner & McDonald 2009; Kumazaki et al. 2011; Meerburg et al. 2012; Ashoorioon & Krause 2006; Ashoorioon et al. 2009), a small number of e-folds (Contaldi et al. 2003; Powell & Kinney 2007; Nicholson & Contaldi 2008), or other exotic inflationary models (Lesgourgues 2000; Feng & Zhang 2003; Mathews et al. 2004; Jain et al. 2009; Romano & Sasaki 2008; Piao et al. 2004; Choudhury et al. 2013; Choudhury & Mazumdar 2014). Therefore, determining the shape of the primordial power spectrum will allow us to evaluate how well these models of the early Universe compare to the observations, rule out some of the proposed models, and thus will give us a better intuition into the conditions of the primordial Universe.

A few probes of the physics of the early Universe include non-Gaussianity, the primordial tensor power spectrum, a cosmic gravitational wave background, and a cosmic neutrino background, none of which have been observed with an acceptable significance. On the other hand, we can observe P(k) through the windows of the cosmic microwave background (CMB) and large scale structure (LSS), which are incredibly important and powerful insights into the early Universe.

The recent Planck mission temperature anisotropy data, combined with the WMAP large-angle polarisation, constrain the scalar spectral index to n_s = 0.9603 ± 0.0073 (Planck Collaboration XXII 2014), which rules out exact scale invariance at over 5σ. In addition, Planck does not find a statistically significant running of the scalar spectral index, obtaining α_s = −0.0134 ± 0.0090. In Planck Collaboration (XXII 2014) an extensive investigation is performed to see whether the primordial power spectrum contains any features. They report that a penalised likelihood approach suggests a feature near the highest wavenumbers probed by Planck at an estimated significance of ~3σ. In addition, a parameterised oscillatory feature does improve the fit to the data by $\mathrm{\Delta }{\mathit{\chi }}_{\mathrm{eff}}^{\mathrm{2}}\mathrm{\approx }\mathrm{10}$, however Bayesian evidence does not prefer these models. On the other hand, high-resolution CMB experiments, such as the South Pole Telescope (SPT)¹, detect a small running of the spectral index; −0.046 < α_s < −0.003 at 95% confidence (Hou et al. 2014). In general, any detections of the running of the spectral index have been small and consistent with zero. Therefore, a highly sensitive algorithm is required to detect these small deviations.

There are generally two approaches to determine the shape of the primordial power spectrum, one by parametrisation and the second by reconstruction. Numerous parametric approaches that search for features with a similar form to those in complex inflationary models have been performed along with a simple binning of P(k) (Bridle et al. 2003; Contaldi et al. 2003; Parkinson et al. 2005; Sinha & Souradeep 2006; Sealfon et al. 2005; Mukherjee & Wang 2005; Bridges et al. 2006; Covi et al. 2006; Hazra et al. 2010; Joy et al. 2009; Verde & Peiris 2008; Paykari & Jaffe 2010; Guo et al. 2011; Goswami & Prasad 2013). Non-parametric methods, which make no assumptions about the model of the early Universe, have also been probed (Hannestad 2001; Wang & Mathews 2002; Matsumiya et al. 2002; Shafieloo & Souradeep 2004; Bridle et al. 2003; Kogo et al. 2004a; Mukherjee & Wang 2003a,b; Hannestad 2004; Kogo et al. 2004b; Tocchini-Valentini et al. 2005; Leach 2006; Shafieloo et al. 2007; Shafieloo & Souradeep 2008; Nagata & Yokoyama 2008, 2009; Nicholson & Contaldi 2009; Nicholson et al. 2010; Hazra et al. 2013). For an extensive review on how to search for features in the primordial power spectrum using a wide range of methods, refer to the following papers and the references therein, which provide a sample on non-parametric reconstruction: deconvolution (Tocchini-Valentini et al. 2006; Ichiki & Nagata 2009; Ichiki et al. 2010), Richardson-Lucy deconvolution (Lucy 1974; Richardson 1972; Hamann et al. 2010; Shafieloo & Souradeep 2008), smoothing splines (Verde & Peiris 2008; Peiris & Verde 2010; Sealfon et al. 2005; Gauthier & Bucher 2012), linear interpolation (Hannestad 2004; Bridle et al. 2003), and Bayesian model selection (Bridges et al. 2009; Vázquez et al. 2012).

Non-parametric methods are hampered by the non-invertibility of the transfer function that descries the transfer from P(k) to CMB (or LSS). Specifically for the CMB power spectrum, the dependence on the transfer function has the form ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathrm{d}\ln \mathit{k}{\mathrm{\Delta }}_{\mathit{\ell }}^{\mathrm{2}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{P}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}$(3)where ℓ is the angular wavenumber that corresponds to an angular scale via ℓ ~ 180°/θ and Δ_ℓ(k) is the angular transfer function of the radiation anisotropies, which holds the cosmological parameters responsible for the evolution of the Universe. As the CMB spectrum is jointly sensitive to the primordial spectrum and the cosmological parameters in the transfer function, there is an induced degeneracy between them. The impact and level of this degeneracy have been investigated in (Paykari & Jaffe 2010). A joint estimation of the cosmological parameters and a free form primordial power spectrum would be prohibitively expensive to perform (as the parameter space can become very large). As a result, a parametric form of the primordial power spectrum is assumed when jointly estimating this spectrum along with the other cosmological parameters. This hides any degeneracies between the cosmological parameters in the transfer function and the form of P(k). Thus it is not clear what the significance of any features found in the reconstructed P(k) should be. One way to break this induced degeneracy is by adding extra information, such as polarisation or LSS data (Hu & Okamoto 2004; Nicholson & Contaldi 2009; Mortonson et al. 2009).

The other hurdle in the estimation of the primordial spectrum is that this continuous spectrum is deconvolved from discrete data C_ℓ. This causes problems if the primordial power spectrum contains features that are smaller or comparable to the gridding in ℓ (Δℓ = 1). This limits our ability to fully recover the primordial power spectrum; in the case of the CMB, even a perfect survey cannot recover the primordial power spectrum completely (Hu & Okamoto 2004).

Here, we propose a new non-parametric method for the reconstruction of the primordial power spectrum from CMB data which is based on the sparsity of the primordial power spectra in a wavelet basis and an appropriate noise modelling of the CMB power spectrum (Paykari et al. 2012).

In Sect. 2 we present the primordial power spectrum reconstruction problem and describe the technique we have developed to perform the reconstruction. Our algorithm is tested on three sets of simulated spectra and applied to WMAP nine-year data in Sect. 3. In Sect. 4 we conclude and give some potential perspectives.

2. Sparse recovery of the primordial power spectrum

2.1. Empirical power spectrum

A CMB experiment, such as Planck, measures the CMB temperature anisotropy Θ(p) in direction p, which is described as T(p) = T_CMB [1 + Θ(p)]. This anisotropy field can be expanded in terms of spherical harmonic functions Y_ℓm as $\mathrm{\Theta }\mathrm{\left(}{p}\mathrm{\right)}\mathrm{=}\sum _{\mathit{\ell }\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\sum _{\mathit{m}\mathrm{=}\mathrm{-}\mathit{\ell }}^{\mathit{\ell }}{\mathit{a}}_{\mathit{\ell m}}{\mathit{Y}}_{\mathit{\ell m}}\mathrm{\left(}{p}\mathrm{\right)}\mathit{,}$(4)with a_ℓm being the spherical harmonic coefficients. The CMB anisotropy Θ(p) is assumed to be Gaussian distributed, which makes the a_ℓm independent and identically distributed (i.i.d.) Gaussian variables with zero mean, ⟨ a_ℓm ⟩ = 0, and variance $\mathrm{⟨}{\mathit{a}}_{\mathit{\ell m}}{\mathit{a}}_{{\mathit{\ell }}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}^{\mathrm{\ast }}\mathrm{⟩}\mathrm{=}{\mathit{\delta }}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}}{\mathit{\delta }}_{\mathit{m}{\mathit{m}}^{\mathrm{\prime }}}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}\mathit{,}$(5)where ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}$ is the CMB temperature angular power spectrum introduced in Eq. (3). However, we only observe a realisation of this underlying power spectrum on our sky, which we can estimate using the empirical power spectrum estimator defined as $\begin{array}{c}\mathrm{􏽢}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}\end{array}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}}\sum _{\mathit{m}\mathrm{=}\mathrm{-}\mathit{\ell }}^{\mathit{\ell }}\mathrm{|}{\mathit{a}}_{\mathit{\ell m}}{\mathrm{|}}^{\mathrm{2}}\mathit{,}$(6)where $\begin{array}{c}{\mathrm{􏽢}}_{}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}\end{array}$ is an unbiased estimator of the true underlying power spectrum; this becomes $\mathrm{⟨}\begin{array}{c}{}_{\mathrm{􏽢}}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}\end{array}\mathrm{⟩}\mathrm{=}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}$ in the case of noiseless CMB data over full sky.

For a given ℓ, the empirical power spectrum follows a χ² distribution with 2ℓ + 1 degrees of freedom, as it is a sum of the squares of independent Gaussian random variables. To account for this variability, we recast the relation between $\begin{array}{c}{\mathrm{􏽢}}_{}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}\end{array}$ and ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}$ as $\begin{array}{c}\mathrm{􏽢}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}\end{array}\mathrm{=}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}{\mathit{Z}}_{\mathit{\ell }}\mathit{,}$(7)where ${\mathit{Z}}_{\mathit{\ell }}\mathrm{=}{\sum }_{\mathit{m}}\mathrm{|}{\mathit{a}}_{\mathit{\ell m}}{\mathrm{|}}^{\mathrm{2}}\mathit{/}\mathit{L}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}$, which is a random variable representing a multiplicative noise distributed according to $\mathit{L}{\mathit{Z}}_{\mathit{\ell }}\mathrm{~}{\mathit{\chi }}_{\mathit{L}}^{\mathrm{2}}\mathit{,} \mathrm{where} \mathit{L}\mathrm{=}\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1.}$(8)In particular, the standard deviation of the empirical power spectrum estimator for a given ℓ is $\sqrt{{}^{\mathrm{(}}\mathrm{2}\mathit{/}{\mathit{L}}^{\mathrm{)}}}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}$.

2.2. Accounting for instrumental noise and partial sky coverage

So far, we have considered that the CMB anisotropy data was available on the full sky which is not possible in practice because of the different Galactic foregrounds. Applying a mask on the sky results in the a modification of the spherical harmonic coefficients of the CMB temperature anisotropy, $\begin{array}{c}\mathrm{˜}\\ {\mathit{a}}_{\mathit{\ell m}}\end{array}\mathrm{=}\mathrm{\int }\mathrm{\Theta }\mathrm{\left(}{p}\mathrm{\right)}\mathit{W}\mathrm{\left(}{p}\mathrm{\right)}{\mathit{Y}}_{\mathit{\ell m}}^{\mathrm{\ast }}\mathrm{\left(}{p}\mathrm{\right)}\mathrm{d}{p}\mathit{,}$(9)where W(p) is the window function applied to the data. The presence of the window function induces correlations between the a_ℓm coefficients at different ℓ and different m and hence Eq. (5) is no longer true.

One can define the pseudo power spectrum $\begin{array}{c}{\mathrm{􏽥}}_{}\\ {\mathit{C}}_{\mathit{\ell }}\end{array}$ as the application of the empirical power spectrum estimator on the spherical harmonic coefficients of the masked sky. When data is contaminated with additive Gaussian stationary noise, the pseudo power spectrum is $\begin{array}{c}\mathrm{􏽥}\\ {\mathit{C}}_{\mathit{\ell }}\end{array}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}}\sum _{\mathit{m}\mathrm{=}\mathrm{-}\mathit{\ell }}^{\mathit{\ell }}\mathrm{|}\begin{array}{c}\mathrm{˜}\\ {\mathit{a}}_{\mathit{\ell m}}\end{array}\mathrm{+}\begin{array}{c}\mathrm{˜}\\ {\mathit{n}}_{\mathit{\ell m}}\end{array}{\mathrm{|}}^{\mathrm{2}}\mathit{,}$(10)where $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{n}}_{\mathit{\ell m}}\end{array}$ are the spherical harmonic coefficients of the masked instrumental noise.

Following the MASTER method from Hivon et al. (2002), the pseudo power spectrum $\begin{array}{c}{\mathrm{􏽥}}_{}\\ {\mathit{C}}_{\mathit{\ell }}\end{array}$ and the empirical power spectrum $\begin{array}{c}{\mathrm{􏽢}}_{}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}\end{array}$ can be related through their ensemble averages, $\mathrm{⟨}\begin{array}{c}\mathrm{􏽥}\\ {\mathit{C}}_{\mathit{\ell }}\end{array}\mathrm{⟩}\mathrm{=}\sum _{{\mathit{\ell }}^{\mathrm{\prime }}}{\mathit{M}}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}}\mathrm{⟨}\begin{array}{c}\mathrm{􏽢}\\ {\mathit{C}}_{{\mathit{\ell }}^{\mathrm{\prime }}}^{\mathrm{th}}\end{array}\mathrm{⟩}\mathrm{+}\mathrm{⟨}\begin{array}{c}\mathrm{􏽥}\\ {\mathit{N}}_{\mathit{\ell }}\end{array}\mathrm{⟩}\mathit{,}$(11)where M_ℓℓ′ describes the mode-mode coupling between modes ℓ and ℓ′ resulting from computing the transform on the masked sky. We note that in this expression $\mathrm{⟨}\begin{array}{c}{}_{\mathrm{􏽢}}\\ {\mathit{C}}_{{\mathit{\ell }}^{\mathrm{\prime }}}^{\mathrm{th}}\end{array}\mathrm{⟩}\mathrm{=}{\mathit{C}}_{{\mathit{\ell }}^{\mathrm{\prime }}}^{\mathrm{th}}$ and we introduce the notations ${\mathit{C}}_{\mathit{\ell }}\mathrm{=}\mathrm{⟨}\begin{array}{c}\mathrm{􏽥}\\ {\mathit{C}}_{\mathit{\ell }}\end{array}\mathrm{⟩} \mathrm{and} {\mathit{N}}_{\mathit{\ell }}\mathrm{=}\mathrm{⟨}\begin{array}{c}\mathrm{􏽥}\\ {\mathit{N}}_{\mathit{\ell }}\end{array}\mathrm{⟩}\mathit{,}$(12)where C_ℓ and N_ℓ refer to the CMB and the noise power spectra of the masked maps, respectively.

We will also work under the approximation that the pseudo power spectrum $\begin{array}{c}{\mathrm{􏽥}}_{}\\ {\mathit{C}}_{\mathit{\ell }}\end{array}$ still follows a χ² distribution with 2ℓ + 1 degrees of freedom and can be modelled as $\begin{array}{ccc}\begin{array}{c}\mathrm{􏽥}\\ {\mathit{C}}_{\mathit{\ell }}\end{array}& \mathrm{=}& {\mathit{C}}_{\mathit{\ell }}{\mathit{Z}}_{\mathit{\ell }}\mathit{,}\\ & \mathrm{=}& \left(\sum _{{\mathit{\ell }}^{\mathrm{\prime }}}{\mathit{M}}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}}{\mathit{C}}_{{\mathit{\ell }}^{\mathrm{\prime }}}^{\mathrm{th}}\mathrm{+}{\mathit{N}}_{\mathit{\ell }}\right){\mathit{Z}}_{\mathit{\ell }}\mathit{,}\end{array}$where Z_ℓ is defined in Eq. (8).

2.3. Formulation of the inverse problem

Here we aim to estimate the primordial power spectrum P_k from the pseudo power spectrum $\begin{array}{c}{\mathrm{􏽥}}_{}\\ {\mathit{C}}_{\mathit{\ell }}\end{array}$ computed on a masked noisy map of the sky.

Equation (13) relates the observables $\begin{array}{c}{\mathrm{􏽥}}_{}\\ {\mathit{C}}_{\mathit{\ell }}\end{array}$ to the theoretical CMB anisotropy power spectrum ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}$, taking into account instrumental noise, sample variance, and masking. The theoretical power spectrum ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}$ is itself related to the primordial power spectrum through the convolution operation defined in Eq. (3). For a finite sampling of the wavenumber k, this convolution can be recast as a matrix operator T acting on the discretely sampled primordial spectrum, now referred to as P_k, ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}\mathrm{\simeq }\sum _{\mathit{k}}{\mathit{T}}_{\mathit{\ell k}}{\mathit{P}}_{\mathit{k}}\mathit{,}$(15)with matrix elements ${\mathit{T}}_{\mathit{\ell k}}\mathrm{=}\mathrm{4}\mathit{\pi }\mathrm{\Delta }\ln \mathit{k}\hspace{0.17em}{\mathrm{\Delta }}_{\mathit{\ell k}}^{\mathrm{2}}$, where Δlnk is the logarithmic k interval for the discrete sampling chosen in the integration of the system of equations. Because of the non-invertibility of the T operator, recovering the primordial power spectrum P_k from the true CMB power spectrum ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}$ constitutes an ill-posed inverse problem. Finally, the complete problem we aim to solve can be condensed in the following form: $\begin{array}{c}\mathrm{􏽥}\\ {\mathit{C}}_{\mathit{\ell }}\end{array}\mathrm{=}\left(\sum _{{\mathit{\ell }}^{\mathrm{\prime }}\mathit{k}}{\mathit{M}}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}}{\mathit{T}}_{{\mathit{\ell }}^{\mathrm{\prime }}\mathit{k}}{\mathit{P}}_{\mathit{k}}\mathrm{+}{\mathit{N}}_{\mathit{\ell }}\right){\mathit{Z}}_{\mathit{\ell }}\mathit{.}$(16)We assume that the masked instrumental noise power spectrum N_ℓ is known for a given experiment. It can be computed from a JackKnife data map or from realistic instrumental noise simulations. Therefore, in the power spectrum of the data $\begin{array}{c}{\mathrm{􏽥}}_{}\\ {\mathit{C}}_{\mathit{\ell }}\end{array}$, only the primordial power spectrum P_k remains unknown. Here we assume that the cosmology is known and hence operator T is known.

The presence of the multiplicative noise Z_ℓ further complicates the ill-posed inverse problem of Eq. (15). We address both the inversion problem and the control of the noise in the framework of sparse recovery. The inversion problem in Eq. (16) can be regularised in a robust way by using the sparse nature of the reconstructed signal as a prior. Furthermore, sparse recovery has already been successfully used in the TOUSI algorithm (Paykari et al. 2012) to handle the multiplicative noise term and denoise the CMB power spectrum with high accuracy from single realisations.

2.4. The TOUSI method

It was shown in Paykari et al. (2012) that the theoretical power spectrum ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}$ can be represented with only a few coefficients (i.e. sparse representation) in a given dictionary (wavelet, DCT, etc.) and that a sparse regularisation allows us to recover the theoretical power spectrum directly from the measured CMB empirical power spectrum $\begin{array}{c}{\mathrm{􏽢}}_{}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}\end{array}$, without having to know the cosmological parameters.

A proper treatment of the non-Gaussian noise on $\begin{array}{c}{\mathrm{􏽢}}_{}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}\end{array}$ was proposed in TOUSI, which is based on the Wahba variance stabilisation transform (VST). After the variance stabilisation is applied, the noise on $\begin{array}{c}{\mathrm{􏽢}}_{}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}\end{array}$ can be treated as an additive Gaussian noise with zero mean and unit variance. The VST operator is defined as $\mathrm{𝒯:}\mathit{x}\mathrm{\in }{R}^{\mathrm{+}}\mathrm{}\mathrm{\to }\frac{\ln \mathit{x}\mathrm{-}{\mathit{\mu }}_{\mathit{L}}}{{\mathit{\sigma }}_{\mathit{L}}}\mathit{,}$(17)where μ_L = ψ₀(L/2) − ln(L/ 2) and ${\mathit{\sigma }}_{\mathit{L}}^{\mathrm{2}}\mathrm{=}{\mathit{\psi }}_{\mathrm{1}}\mathrm{(}\mathit{L}\mathit{/}\mathrm{2}\mathrm{)}$, where ψ_m is the polygamma function ${\mathit{\psi }}_{\mathit{m}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\frac{{\mathrm{d}}^{\mathit{m}\mathrm{+}\mathrm{1}}}{\mathrm{d}{\mathit{t}}^{\mathit{m}\mathrm{+}\mathrm{1}}}\ln \mathrm{\Gamma }\mathrm{\left(}\mathit{t}\mathrm{\right)}$. We denote ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{s}}$ as the stabilised empirical power spectrum after applying the VST and get ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{s}}\mathrm{=}\mathrm{𝒯}\mathrm{\left(}\begin{array}{c}\mathrm{􏽢}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}\end{array}\mathrm{\right)}\mathrm{=}\frac{\ln {\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}}{{\mathit{\sigma }}_{\mathit{L}}}\mathrm{+}{\mathit{ϵ}}_{\mathit{\ell }}\mathit{,}$(18)where . We define the inverse operator of as $\mathrm{ℛ:}\mathit{x}\mathrm{\in }R\mathrm{}\mathrm{\to }\exp \mathrm{\left(}{\mathit{\sigma }}_{\mathit{L}}\mathit{x}\mathrm{\right)}\mathit{.}$(19)Having X_ℓ as the unknown power spectrum to be recovered, the TOUSI method consists of minimising the constrained optimisation problem $\underset{{\mathit{X}}_{\mathit{\ell }}}{\min }\mathrm{\parallel }{{\Phi }}^{\mathrm{t}}{\mathit{X}}_{\mathit{\ell }}{\mathrm{\parallel }}_{\mathrm{1}} \mathrm{s}\mathit{.}\mathrm{t}\mathit{.} \{\begin{array}{c}\\ \\ \end{array}$(20)where ${\mathit{Y}}_{\mathit{\ell }}\mathrm{=}{\mathit{X}}_{\mathit{\ell }}\mathrm{+}{\mathit{N}}_{\mathit{\ell }}^{\mathrm{th}}$; ⊙ stands for the Hadamard product (i.e. entry-wise multiplication) of two vectors; and Φ is the chosen dictionary. Vector S provides a set of active coefficients (not due to noise), where S_i = 1 if the ith coefficient is above the noise level (i.e. significant) and 0 otherwise. This minimisation is performed iteratively, $\begin{array}{ccc}\begin{array}{c}\mathrm{􏽥}\\ {\mathit{X}}_{\mathit{\ell }}\end{array}& \mathrm{=}& \mathrm{ℛ}\mathrm{(}\mathrm{𝒯}\mathrm{(}{{\mathit{Y}}_{\mathit{\ell }}^{\mathit{n}}}^{\mathrm{)}}\mathrm{+}{\Phi }\mathit{S}\mathrm{\odot }{\mathrm{(}{{\Phi }}^{\mathrm{t}}{\mathrm{(}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{s}}\mathrm{-}\mathrm{𝒯}{\mathrm{(}{{\mathit{Y}}_{\mathit{\ell }}^{\mathit{n}}}^{\mathrm{)}}}^{\mathrm{)}}}^{\mathrm{)}}}^{\mathrm{)}}\mathrm{-}{\mathit{N}}_{\mathit{\ell }}^{\mathrm{th}}\mathit{,}\\ {\mathit{X}}_{\mathit{\ell }}^{\mathit{n}\mathrm{+}\mathrm{1}}& \mathrm{=}& {\mathrm{𝒫}}_{\mathrm{+}}\mathrm{(}{\Phi }{\mathrm{ST}}_{{\mathit{\lambda }}_{\mathit{n}}}\mathrm{(}{{\Phi }}^{\mathrm{t}}\begin{array}{c}\mathrm{􏽥}\\ {\mathit{X}}_{\mathit{\ell }}\end{array}{\mathrm{)}}^{\mathrm{)}}\mathit{,}\end{array}$(21)where n is the iteration number, is a positivity constraint. The soft thresholding operator ST_λn has an iteration dependent threshold level λ_n and is defined as $\mathrm{\forall }{x}\mathrm{\in }{R}^{\mathit{n}}\mathit{,}{S}{{T}}_{\mathit{\lambda }}\mathrm{(}{x}{\mathrm{)}}_{\mathit{i}}\mathrm{=}\mathrm{sgn}\mathrm{(}{\mathit{x}}_{\mathit{i}}\mathrm{\left)}\mathrm{\right(}\mathrm{\left|}{\mathit{x}}_{\mathit{i}}\mathrm{\right|}\mathrm{-}\mathit{\lambda }{\mathrm{)}}_{\mathrm{+}}\mathit{.}$(22)Full details of the TOUSI algorithm can be found in Paykari et al. (2012).

2.5. Pk sparse recovery formulation

The problem of reconstructing the primordial power spectrum is stated in Eq. (16). Solving this problem has three inherent difficulties: 1) the singularity of the convolution operator T_ℓk, which makes the inverse problem ill-posed even in the absence of noise; 2) the multiplicative noise on the power spectrum; and 3) the mask applied to the maps, inducing correlations on the power spectrum.

To address the inverse problem, we adopt the sparse regularisation framework. If the signal to recover, P_k in our case, can be sparsely represented in an adapted dictionary Φ, then this problem, known as the basis pursuit denoising BPDN, can be recast as an optimisation problem. In the case of the inverse problem stated in Eq. (16), the optimisation problem can be formulated as $\underset{\mathit{X}}{\min }\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\parallel }{\mathit{C}}_{\mathit{\ell }}\mathrm{-}\mathrm{(}{MT}\mathit{X}\mathrm{+}{\mathit{N}}_{\mathit{\ell }}\mathrm{)}{\mathrm{\parallel }}_{\mathrm{2}}^{\mathrm{2}}\mathrm{+}\mathit{\lambda }\mathrm{\parallel }{{\Phi }}^{\mathrm{t}}\mathit{X}{\mathrm{\parallel }}_{\mathrm{0}}\mathit{,}$(23)where X is the reconstructed estimate for the primordial power spectrum P_k. The first term in Eq. (23) imposes a ℓ₂ fidelity constraint to the data while the second term promotes the sparsity of the solution in dictionary Φ. The parameter λ tunes the sparsity constraint.

One can notice that in Eq. (23), only the ensemble mean of the pseudo power spectrum C_ℓ appears (which is unknown) and not the actual measurements $\begin{array}{c}{\mathrm{􏽥}}_{}\\ {\mathit{C}}_{\mathit{\ell }}\end{array}$. This is linked to the second difficulty; the measurements are contaminated with a multiplicative noise which cannot be handled with the formulation of Eq. (23). This formalism holds for measurements contaminated with additive Gaussian noise which is not the case of the $\begin{array}{c}{\mathrm{􏽥}}_{}\\ {\mathit{C}}_{\mathit{\ell }}\end{array}$. To overcome this issue, we use the variance stabilisation introduced in the TOUSI algorithm.

We let R_ℓ(X) be the residual between C_ℓ and the reconstructed CMB power spectrum given a primordial power spectrum X, C_ℓ(X) = (MTX + N_ℓ): ${\mathit{R}}_{\mathit{\ell }}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathrm{=}{\mathit{C}}_{\mathit{\ell }}\mathrm{-}{\mathit{C}}_{\mathit{\ell }}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathit{.}$(24)We note that R_ℓ(X) is the data fidelity term in Eq. (23). Since C_ℓ is unknown, so is R_ℓ(X), but we can estimate it from the data $\begin{array}{c}{\mathrm{􏽥}}_{}\\ {\mathit{C}}_{\mathit{\ell }}\end{array}$. We consider the difference: $\begin{array}{ccc}\mathrm{𝒯}\mathrm{\left(}\begin{array}{c}\mathrm{􏽥}\\ {\mathit{C}}_{\mathit{\ell }}\end{array}\mathrm{\right)}\mathrm{-}\frac{\ln \mathrm{\left(}{\mathit{C}}_{\mathit{\ell }}\mathrm{\right(}\mathit{X}\mathrm{\left)}\mathrm{\right)}}{{\mathit{\sigma }}_{\mathit{L}}}& \mathrm{=}& \frac{\ln \mathrm{\left(}{\mathit{C}}_{\mathit{\ell }}\mathrm{\right)}\mathrm{-}\ln \mathrm{\left(}{\mathit{C}}_{\mathit{\ell }}\mathrm{\right(}\mathit{X}\mathrm{\left)}\mathrm{\right)}}{{\mathit{\sigma }}_{\mathit{L}}}\mathrm{+}{\mathit{ϵ}}_{\mathit{\ell }}\mathit{,}\\ & \mathrm{=}& \frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathit{L}}}\ln \left(\frac{{\mathit{C}}_{\mathit{\ell }}}{{\mathit{C}}_{\mathit{\ell }}\mathrm{\left(}\mathit{X}\mathrm{\right)}}\right)\mathrm{+}{\mathit{ϵ}}_{\mathit{\ell }}\mathit{,}\\ & \mathrm{=}& \frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathit{L}}}\ln \left(\mathrm{1}\mathrm{+}\frac{{\mathit{R}}_{\mathit{\ell }}\mathrm{\left(}\mathit{X}\mathrm{\right)}}{{\mathit{C}}_{\mathit{\ell }}\mathrm{\left(}\mathit{X}\mathrm{\right)}}\right)\mathrm{+}{\mathit{ϵ}}_{\mathit{\ell }}\mathit{,}\end{array}$where ϵ_ℓ is the Gaussian noise with zero mean introduced in Eq. (18). Assuming that the residual R_ℓ(X) is small compared to C_ℓ(X), one can linearise the above equation, to a good approximation, as $\mathrm{𝒯}\mathrm{\left(}\begin{array}{c}\mathrm{􏽥}\\ {\mathit{C}}_{\mathit{\ell }}\end{array}\mathrm{\right)}\mathrm{-}\frac{\ln \mathrm{\left(}{\mathit{C}}_{\mathit{\ell }}\mathrm{\right(}\mathit{X}\mathrm{\left)}\mathrm{\right)}}{{\mathit{\sigma }}_{\mathit{L}}}\mathrm{\simeq }\frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathit{L}}{\mathit{C}}_{\mathit{\ell }}\mathrm{\left(}\mathit{X}\mathrm{\right)}}{\mathit{R}}_{\mathit{\ell }}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathrm{+}{\mathit{ϵ}}_{\mathit{\ell }}\mathit{,}$(28)and ${\mathit{R}}_{\mathit{\ell }}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathrm{\simeq }{\mathit{C}}_{\mathit{\ell }}\mathrm{\left(}\mathit{X}\mathrm{\right)}{\mathit{\sigma }}_{\mathit{L}}\left(\mathrm{𝒯}\mathrm{\left(}\begin{array}{c}\mathrm{􏽥}\\ {\mathit{C}}_{\mathit{\ell }}\end{array}\mathrm{\right)}\mathrm{-}\frac{\ln \mathrm{\left(}{\mathit{C}}_{\mathit{\ell }}\mathrm{\right(}\mathit{X}\mathrm{\left)}\mathrm{\right)}}{{\mathit{\sigma }}_{\mathit{L}}}\right)\mathrm{-}{\mathit{C}}_{\mathit{\ell }}\mathrm{\left(}\mathit{X}\mathrm{\right)}{\mathit{\sigma }}_{\mathit{L}}{\mathit{ϵ}}_{\mathit{\ell }}\mathit{.}$(29)In this expression, the variance of the noise, i.e. the second term in the above equation, depends on the current estimate C_ℓ(X). As we need to estimate the variance of the noise propagated to the wavelet coefficients using Monte Carlo simulations, it would be too expensive to estimate this every time C_ℓ(X) changes. Therefore, we opted for an additional approximation and replace the term C_ℓ(X)σ_L by C_ℓ(X⁰)σ_L, where X⁰ is now a fixed fiducial power spectrum which can be the initial guess of the solution. We can now introduce the estimator ${\overline{\mathit{R}}}_{\mathit{\ell }}\mathrm{\left(}\mathit{X}\mathrm{\right)}$ for R_ℓ(X) defined as ${\overline{\mathit{R}}}_{\mathit{\ell }}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathrm{\equiv }{\mathit{C}}_{\mathit{\ell }}\mathrm{\left(}{\mathit{X}}^{\mathrm{0}}\mathrm{\right)}{\mathit{\sigma }}_{\mathit{L}}\left(\mathrm{𝒯}\mathrm{\left(}\begin{array}{c}\mathrm{􏽥}\\ {\mathit{C}}_{\mathit{\ell }}\end{array}\mathrm{\right)}\mathrm{-}\frac{\ln \mathrm{\left(}{\mathit{C}}_{\mathit{\ell }}\mathrm{\right(}\mathit{X}\mathrm{\left)}\mathrm{\right)}}{{\mathit{\sigma }}_{\mathit{L}}}\right)\mathit{,}$(30)which leads to ${\overline{\mathit{R}}}_{\mathit{\ell }}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathrm{\simeq }\frac{{\mathit{C}}_{\mathit{\ell }}\mathrm{\left(}{\mathit{X}}^{\mathrm{0}}\mathrm{\right)}}{{\mathit{C}}_{\mathit{\ell }}\mathrm{\left(}\mathit{X}\mathrm{\right)}}{\mathit{R}}_{\mathit{\ell }}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathrm{+}{\mathit{C}}_{\mathit{\ell }}\mathrm{\left(}{\mathit{X}}^{\mathrm{0}}\mathrm{\right)}{\mathit{\sigma }}_{\mathit{L}}{\mathit{ϵ}}_{\mathit{\ell }}\mathit{.}$(31)Unless C_ℓ(X⁰) = C_ℓ(X) in the first term, this estimator yields a biased estimate of the amplitude of R_ℓ(X). However, it still verifies ${\overline{\mathit{R}}}_{\mathit{\ell }}\mathrm{\left(}{\mathit{P}}_{\mathit{k}}^{\mathrm{th}}\mathrm{\right)}\mathrm{=}\mathrm{0}$ and unless the estimated solution X deviates significantly from X⁰, the ratio C_ℓ(X⁰) /C_ℓ(X) remains limited to within a few percents. Furthermore, the fiducial power spectrum X⁰ can be reset several times to the current estimated X as the algorithm converges towards a solution, thereby removing any potential multiplicative bias on the residuals once the algorithm has converged. On the other hand, the noise on the estimator ${\overline{\mathit{R}}}_{\mathit{\ell }}\mathrm{\left(}\mathit{X}\mathrm{\right)}$ now has a fixed variance independent of the current estimate of the solution X. Replacing this estimator in the data fidelity term of Eq. (23) eliminates the unknown true anisotropy power spectrum from the data fidelity term.

Furthermore, we modify the sparsity constraint by applying a weight for each wavelet coefficient, thus turning the parameter λ in Eq. (23) into Kλ_i, where i is the coefficient index in the wavelet domain. In Sect. 2.7, a specific choice of the λ_i will allow us to use a single regularisation parameter K to handle the non-stationary and correlated noise on the estimator ${\overline{\mathit{R}}}_{\mathit{\ell }}$ in a way that translates into a significance level threshold for the detection of features. The optimisation problem solved by PRISM can now be formulated as $\underset{\mathit{X}}{\min }\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\parallel }\frac{\mathrm{1}}{{\mathit{C}}_{\mathit{\ell }}\mathrm{\left(}{\mathit{X}}^{\mathrm{0}}\mathrm{\right)}{\mathit{\sigma }}_{\mathit{L}}}{\overline{\mathit{R}}}_{\mathit{\ell }}\mathrm{\left(}\mathit{X}\mathrm{\right)}{\mathrm{\parallel }}_{\mathrm{2}}^{\mathrm{2}}\mathrm{+}\mathit{K}\sum _{\mathit{i}}{\mathit{\lambda }}_{\mathit{i}}\mathrm{\parallel }\mathrm{[}{{\Phi }}^{\mathrm{t}}\mathit{X}{\mathrm{]}}_{\mathit{i}}{\mathrm{\parallel }}_{\mathrm{0}}\mathit{,}$(32)where the pre-factor 1 /C_ℓ(X⁰)σ_L weights the ℓ₂ data fidelity term according to the variance of the noise on the estimator ${\overline{\mathit{R}}}_{\mathit{\ell }}$.

2.6. The PRISM algorithm

The ℓ₀ optimisation problem stated in Eq. (32)cannot be solved directly. However, the solution can be estimated by solving a sequence of relaxed problems using the re-weighted ℓ₁ minimisation technique Candes et al. (2008). This technique amounts to solving a sequence of weighted ℓ₁ problems of the form $\underset{\mathit{X}}{\min }\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\parallel }\frac{\mathrm{1}}{{\mathit{C}}_{\mathit{\ell }}\mathrm{\left(}{\mathit{X}}^{\mathrm{0}}\mathrm{\right)}{\mathit{\sigma }}_{\mathit{L}}}{\overline{\mathit{R}}}_{\mathit{\ell }}\mathrm{\left(}\mathit{X}\mathrm{\right)}{\mathrm{\parallel }}_{\mathrm{2}}^{\mathrm{2}}\mathrm{+}\mathit{K}\sum _{\mathit{i}}{\mathit{\lambda }}_{\mathit{i}}\mathrm{\left|}\mathrm{\right[}{W}{{\Phi }}^{\mathrm{t}}\mathit{X}{\mathrm{]}}_{\mathit{i}}\mathrm{|}\mathit{,}$(33)where W is a diagonal matrix applying a different weight for each wavelet coefficient. This relaxed problem is now tractable and the solution of the original problem (32)can be approximated using the iterative algorithm presented in Candes et al. (2008) to perform the reweighted analysis-based ℓ₁ recovery:

• 1.

  Set j = 0, for each element of the weighting matrix W set ${\mathit{w}}_{\mathit{i}}^{\mathit{j}}\mathrm{=}\mathrm{1}$. Set the first guess X⁰ by fitting a pure scale-invariant primordial power spectrum to the data $\begin{array}{c}{\mathrm{􏽥}}_{}\\ {\mathit{C}}_{\mathit{\ell }}\end{array}$.

• 2.

  Solve the weighted ℓ₁ problem (33)yielding a solution X^j.

• 3.

  Compute ${\mathit{\alpha }}_{\mathit{i}}^{\mathit{j}}\mathrm{=}{\Phi }{\mathit{X}}^{\mathit{j}}$ and update the weights according to: ${\mathit{w}}_{\mathit{i}}^{\mathit{j}\mathrm{+}\mathrm{1}}\mathrm{=}\{\begin{array}{c}\\ \frac{\mathrm{1}}{\mathrm{\left|}{\mathit{\alpha }}_{\mathit{i}}^{\mathit{j}}\mathrm{\right|}\mathit{/}\mathit{K}{\mathit{\lambda }}_{\mathit{i}}}& \mathrm{if}\mathrm{\left|}{\mathit{\alpha }}_{\mathit{i}}^{\mathit{j}}\mathrm{\right|}\mathrm{\ge }\mathit{K}{\mathit{\lambda }}_{\mathit{i}}\\ \mathrm{1}& \mathrm{if}\mathrm{\left|}{\mathit{\alpha }}_{\mathit{i}}^{\mathit{j}}\mathrm{\right|}\mathit{<}\mathit{K}{\mathit{\lambda }}_{\mathit{i}}\mathit{,}\end{array}$(34)where λ_i is the standard deviation propagated to the wavelet coefficients (see Sect. 2.7) and K is a given significance level.

• 4.

  Terminate on convergence or when reaching the maximum number of iterations, otherwise go to step 2.

In practice, we find that three iterations of this procedure are enough to reach satisfying convergence and de-biasing of our results and we see no further improvements by performing additional re-weightings.

To solve the relaxed problem (33)given a weighting matrix W, the popular iterative soft-thresholding algorithm (ISTA) can be used. This proximal forward-backward iterative scheme relies on the iteration $\begin{array}{ccc}\begin{array}{c}\mathrm{􏽥}\\ {\mathit{X}}^{\mathit{n}\mathrm{+}\mathrm{1}}\end{array}& \mathrm{=}& {\mathit{X}}^{\mathit{n}}\mathrm{+}\mathit{\mu }{{T}}^{\mathrm{t}}{{M}}^{\mathrm{t}}\frac{\mathrm{1}}{\mathrm{\left(}{\mathit{C}}_{\mathit{\ell }}\mathrm{\right(}{\mathit{X}}^{\mathrm{0}}\mathrm{)}{\mathit{\sigma }}_{\mathit{L}}{\mathrm{)}}^{\mathrm{2}}}{\overline{\mathit{R}}}_{\mathit{\ell }}\mathrm{\left(}{\mathit{X}}^{\mathit{n}}\mathrm{\right)}\mathit{,}\\ {\mathit{X}}^{\mathit{n}\mathrm{+}\mathrm{1}}& \mathrm{=}& {\mathrm{prox}}_{\mathit{K\mu }\mathrm{\parallel }\mathit{\lambda }\mathrm{\odot }\mathit{W}{\mathrm{\Phi }}^{\mathrm{t}}\mathrm{·}{\mathrm{\parallel }}_{\mathrm{1}}}\mathrm{(}{\begin{array}{c}\mathrm{􏽥}\\ {\mathit{X}}^{\mathit{n}\mathrm{+}\mathrm{1}}\end{array}}^{\mathrm{)}}\mathit{,}\end{array}$where μ is an adapted step size and prox_{Kμ ∥ λ⊙WΦ^t· ∥ 1} is the proximal operator corresponding to the sparsity constraint. The gradient descent step μ has to verify $\mathrm{0}\mathit{<}\mathit{\mu }\mathrm{\le }\frac{\mathrm{2}}{\mathrm{\parallel }{{T}}^{\mathrm{t}}{{M}}^{\mathrm{t}}\mathrm{\left(}{\mathit{C}}_{\mathit{\ell }}\mathrm{\right(}{\mathit{X}}^{\mathrm{0}}\mathrm{)}{\mathit{\sigma }}_{\mathit{L}}{\mathrm{)}}^{-2}{MT}\mathrm{\parallel }}\mathit{,}$(37)where ∥·∥ is the spectral norm of the operator.

In the absence of a closed-form expression for the proximal operator, its value can be estimated by solving a nested optimisation problem: $\{\begin{array}{c}\\ \mathit{û}\mathrm{=}{\arg \hspace{0.17em}\min }_{\mathrm{\left|}{\mathit{u}}_{\mathit{i}}\mathrm{\right|}\mathrm{\le }\mathit{K\mu }{\mathit{\lambda }}_{\mathit{i}}{\mathit{w}}_{\mathit{i}}}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\parallel }{\Phi }\mathit{u}\mathrm{-}\mathit{x}{\mathrm{\parallel }}_{\mathrm{2}}^{\mathrm{2}}\\ {\mathrm{prox}}_{\mathit{K\mu }\mathrm{\parallel }\mathit{\lambda }\mathrm{\odot }\mathit{W}{{\Phi }}^{\mathrm{t}}\mathrm{·}{\mathrm{\parallel }}_{\mathrm{1}}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\mathit{x}\mathrm{-}{\Phi }\mathit{û}\mathit{.}\end{array}$(38)We solve this optimisation problem at each iteration of the algorithm, using the fast iterative soft-thresholding algorithm (FISTA) Beck & Teboulle (2009), a fast variant of ISTA. The details of the algorithm solving this weighted problem are provided in Algorithm 1.

2.7. Choice of wavelet dictionary and regularisation parameter

As mentioned in the previous section, the regularisation parameter K can be set according to a desired significance level. In Eq. (38), it can be seen that the wavelet coefficients u_i are constrained within a weighted ℓ₁ ball and correspond to the non-significant part of the signal. In order to place the radius of this ℓ₁ ball according to the expected level of noise for each wavelet coefficient, we propagate the noise on the estimator ${\overline{\mathit{R}}}_{\mathit{\ell }}$ from Eq. (35)through the operator ΦT^tM^t(C_ℓ(X⁰)σ_L)^-2 and estimate its standard deviation at each pixel and each wavelet scale. In practice, we estimate this noise level using Monte Carlo simulations of the noise on ${\overline{\mathit{R}}}_{\mathit{\ell }}$. We set each λ_i to the resulting standard deviation for each wavelet coefficient. As a result, coefficients below Kλ_i will be considered as part of the noise and one only need to set a global parameter K to tune the sparsity constraint according to the noise level. In the following section, we have chosen to set this parameter at K = 5, thus robustly removing noise.

The choice of wavelet Φ will have an impact on the performance of the algorithm. In the following study, we use bi-orthogonal Battle-Lemarié wavelets of order 1 with 9 dyadic wavelet scales. This choice of wavelet is generic and not specifically tuned to a type of primordial power spectrum. More physically motivated dictionaries could be used to reconstruct a specific type of feature predicted by a given theory.

3. Results

3.1. Numerical simulations

To assess the performance of our non-linear Algorithm 1 we perform a series of reconstructions for three different types of primordial power spectra: a near scale-invariant spectrum with n_s = 0.972 (Hinshaw et al. 2013), a spectrum with a small running of the spectral index with n_s = 0.972 and α_s = −0.017 (Hou et al. 2014), and a spectrum with n_s = 0.972 with a compensated feature around k = 0.03 Mpc^-1. The first two simple models are the most favoured by the current data and the spectrum with the feature (investigated in other works, see Nicholson & Contaldi 2009) is only used to demonstrate the ability of the algorithm to detect and reconstruct isolated features. In all cases, the cosmological parameters responsible for the evolution of the Universe in the radiation transfer function are kept the same and according to the WMAP nine-year parameters (Hinshaw et al. 2013), Ω_bh² = 0.02264, Ω_ch² = 0.1138, Ω_Λ = 0.721, and τ = 0.089.

For a thorough comparison of our simulations to the WMAP nine-year data we perform the Monte Carlo simulations at the level of the five WMAP frequency channels, taking into account the propagation of the instrumental noise through the component separation and masking steps. For each of the three test primordial spectra we produce a set of 2000 pseudo power spectra $\begin{array}{c}{\mathrm{􏽥}}_{}\\ {\mathit{C}}_{\mathit{\ell }}\end{array}$ by processing the simulated channels through the LGMCA component separation pipeline (Bobin et al. 2013) before computing the empirical power spectrum of the masked maps. In detail, the simulations are produced using the following steps.

• Frequency channels: we simulate CMB maps at the five WMAPchannels at frequencies 23, 33, 41, 61, and 94 GHz. The frequencydependant beams are perfectly isotropic PSFs and their profileshave been obtained as the mean value of the beam transferfunctions at each frequency as provided by the WMAPconsortium (nine year version).

• Instrumental noise: noise maps for each channel have been generated as Gaussian realisations of pixel variance maps obtained by combining the nine one-year full-resolution hit maps as provided by the WMAP consortium.

• Cosmic microwave background: Gaussian realisations of the CMB are computed from the three power spectra ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}$, which were obtained by applying the radiation transfer function T to each of the three test primordial power spectra. The transfer function is computed using CLASS² (Blas et al. 2011) according to the best-fit WMAP nine-year cosmology. The CMB signal for each channel is then obtained by applying the corresponding beam to the simulated CMB map as well as the HEALPix window for nside of 1024.

• LGMCA Component Separation: full sky 15 arcmin resolution maps are obtained by applying LGMCA, with the precomputed set of parameters (Bobin et al. 2013), to the five simulated channels for CMB and noise. Noisy full sky maps are obtained by adding the resulting signal and noise maps.

• Masking: final maps are obtained by applying the WMAP mask kq85 mask with f_sky = 0.75.

The pseudo power spectra are obtained by applying the empirical power spectrum estimator to the simulated maps. The noise power spectrum N_ℓ is estimated by averaging the 2000 pseudo spectra of masked noise maps. Figure 1 shows an example of a masked noisy CMB map obtained from our simulation process. Figure 2 shows the pseudo power spectra for the three test primordial spectra as well as the instrumental noise power spectrum estimated from the simulations. The light blue crosses show one realisation of the pseudo power spectrum for the near scale-invariant primordial power spectrum and the pink crosses show the one with a small running. As can be seen, the three different CMB spectra lie well within each other’s noise band and on large and small scales they become almost indistinguishable. Hence to accurately reconstruct the three underlying primordial power spectra from these CMB spectra, a very good handle on both the instrumental noise and the sample variance is required.

Fig. 1 A simulated noisy CMB map at 15 arcmin resolution obtained from LGMCA and masked with the WMAP kq85 mask. The noise level corresponds to the WMAP nine-year data. This map was generated from a CMB power spectrum for a primordial spectrum with ns = 0.972 and αs = 0.

Fig. 2 CMB pseudo power spectra for the three types of primordial power spectra. The solid blue line shows the pseudo spectrum based on a primordial spectrum with ns = 0.972 and αs = 0. The light blue crosses show one simulation of this spectrum, computed from the map in Fig. 1. The red line shows the pseudo spectrum for a primordial spectrum with ns = 0.972 and αs = − 0.017 and the orange line corresponds to a power spectrum with a localised feature at k = 0.03 Mpc-1. These spectra include the effects of the mask, the 15 arcmin beam, the HEALPix window for nside of 1024, and the instrumental noise power spectrum, which is shown by a solid black line.

3.2. Reconstructions of primordial power spectra

To apply PRISM to the simulated data, we build a transfer function T′ adapted to the simulations so that it includes the effects of the 15 arcmin beam from LGMCA and the HEALPix window of nside = 1024. Using the same radiation transfer function T as computed for the simulations, the resulting transfer matrix T′ can be written as ${{T}}^{\mathrm{\prime }}\mathrm{=}{BHTQ}\mathit{,}$(39)where ${B}\mathrm{=}\mathrm{diag}\mathrm{\left(}{\mathit{b}}_{\mathit{\ell }}^{\mathrm{2}}\mathrm{\right)}$ and ${H}\mathrm{=}\mathrm{diag}\mathrm{\left(}{\mathit{h}}_{\mathit{\ell }}^{\mathrm{2}}\mathrm{\right)}$ with b_ℓ and h_ℓ being the beam and the HEALPix window, respectively. Operator Q performs a linear interpolation from a logarithmic scale using 838 points to the linear sampling in k of the CLASS transfer function T in the range k ~ 10^-4−0.15 Mpc^-1. We also compute the MASTER coupling matrix M_kq85 corresponding to the kq85 high-resolution temperature analysis mask used in the simulations.

We now have all the ingredients necessary in our algorithm: M_kq85, T′, and Φ, which we use to construct our algorithm and apply it to the 3 × 2000 simulated pseudo power spectra. We use the same set of hyper parameters in PRISM for three types of primordial spectra: a Kσ significance level for the sparsity constraint with K = 5, three reweightings, and N_max = 400 iterations per reweighting.

Fig. 3 Reconstructions for the primordial power spectra and their corresponding CMB pseudo spectra. In blue we show the 2000 reconstructed spectra with ns = 0.972 and αs = 0 and in cyan the reconstruction for ns = 0.972 and αs = −0.017. In both cases the mean of the reconstructions is shown in orange and the fiducial input spectrum is shown in red. As can be seen, for k> 0.015Mpc-1 PRISM can reconstruct the primordial power spectra with such accuracy that the two are easily distinguishable, despite their very similar forms in Cℓ space. The shaded regions in panel b) correspond to the 1σ sample (cosmic) variance, which demonstrates the similarity of the two types of CMB spectra. The quality of the reconstruction can also be seen in the reconstructed angular power spectra which are extremely close to the theory and well within the 1σ sample variance intervals.

In Fig. 3a we show the reconstructed primordial spectra in the range k ~ 0.001−0.10Mpc^-1. The blue lines show the 2000 reconstructed spectra for the spectrum with n_s = 0.972 and α_s = 0.0 and the cyan lines show the reconstructions for the spectrum with n_s = 0.972 and α_s = −0.017. In each case, the orange line is the mean of the reconstructions and the red line is the fiducial one.

The reconstruction of the primordial power spectrum is limited by different effects on different scales. On very large scales, there are fundamental physical limitations placed on the recovery of the primordial power spectrum by both the cosmic variance and the more severe geometrical projection of the modes. The physical limitations in the radiation transfer function places an inherent limitation at large scales meaning the primordial power spectrum cannot be fully recovered on these scales, even in a perfect CMB measurement. On the other hand, on small scales we are limited by the instrumental noise. This leaves us with a window through which we can recover the primordial power spectrum with a good accuracy. Nevertheless, as can be seen, for k> 0.015Mpc^-1 the PRISM algorithm can reconstruct the primordial power spectrum to a great accuracy and easily distinguishes between the two types of spectra.

Figure 3b shows the 2000 CMB spectra obtained from the reconstructed primordial power spectra of each type. The blue lines show the CMB power spectra obtained from the near scale-invariant primordial spectra and the cyan lines show the ones for the primordial spectrum with a running. In each case, the orange line shows the mean of the reconstructions and the red line shows the fiducial spectrum. Comparing these CMB spectra to the input simulated ones, shown in Fig. 2, shows the great performance of the PRISM algorithm.

Figure 4 shows the performance of PRISM in reconstructing a localised feature in the primordial power spectrum. The green lines show the 2000 individual reconstructions, the solid orange line shows the mean of the reconstructions, and the fiducial spectrum is shown in red. As can be seen, both the position and the amplitude of the feature can be recovered with great accuracy.

3.3. Reconstruction from WMAP nine-year CMB spectrum

In the WMAP nine-year analysis (Hinshaw et al. 2013), the cosmological parameters in the radiation transfer function are fitted along with n_s and A_s, hence a power law form for the primordial power spectrum is assumed. This means the transfer function computed using these best-fit parameters will always allow a power-law primordial power spectrum to fit the observed data. However, reconstructing a free form primordial power spectrum from the data, assuming the fiducial transfer function, allows us to test this null hypothesis by looking for significant deviations between the reconstructed spectrum from data and the simulations.

The WMAP nine-year data is processed using LGMCA as described in Bobin et al. (2013), which is the same pipeline used to produce the simulations. As mentioned previously, a good handle on the noise power spectrum is critical in order to yield an unbiased reconstruction of the primordial power spectrum. We estimate the noise power spectrum from the WMAP nine-year data by subtracting the cross-power spectrum from the auto-power spectrum and applying a denoising, using the TOUSI algorithm. To account for the effect of point sources, which were not accounted for in the simulations, we add an estimate of the point sources power spectrum, computed from 100 simulations, to the estimated noise power spectrum. Figure 5b shows the pseudo-power spectrum computed from the LGMCA WMAP nine-year map (blue crosses) and the estimated instrumental noise power spectrum (solid black line). We note that in theory, the noise power spectrum could be computed from simulations. However, after comparing our estimated noise power spectrum from the 2000 simulations to the actual noise power spectrum in the WMAP nine-year data we found a small bias that we could not account for in the simulations. Hence we opted to use the data itself to estimate the noise power spectrum.

We apply PRISM, with the same hyper parameters as in the simulations, to the WMAP nine-year LGMCA CMB pseudo power spectrum. The reconstructed primordial power spectrum is shown in red in Fig. 5a. In this figure, we overlay the 1σ interval around the mean of reconstructed primordial near scale-invariant spectrum, obtained from the simulations. The best-fit power-law power spectrum from WMAP nine-year data with n_s = 0.972 and α_s = 0 is shown in yellow, while the best-fit power spectrum with a running from WMAP nine-year data with n_s = 1.009 and α_s = − 0.019 is shown in cyan (Hinshaw et al. 2013). As can be seen, the reconstructed power spectrum from the data does not exhibit a significant deviation from the best-fit near scale-invariant spectrum. The small departure from the 1σ interval at small scales is not significant, especially since our simulations did not thoroughly take into account additional effects such as a beam uncertainty and point sources. To conclude, we find no significant departure from the WMAP nine-year best-fit near scale-invariant spectrum.

Fig. 4 Reconstruction of the primordial power spectrum with ns = 0.972, αs = 0.0, and an additional feature around k = 0.03 Mpc-1 shown in green. The 2000 reconstructions are superimposed with their mean shown in orange. The fiducial input spectrum is shown in red. As can be seen, PRISM is able to recover both the position and the amplitude of the feature with great accuracy.

Fig. 5 Reconstruction of the primordial power spectrum from the LGMCA WMAP nine-year data and its corresponding pseudo spectrum shown in red. For comparison, we also show the mean of the reconstruction for ns = 0.972 and αs = 0 with a solid dark blue line with the 1σ interval around the mean shown as a shaded blue region. The WMAP nine-year fiducial primordial power spectrum with ns = 0.972 and αs = 0 is shown in yellow and in cyan we show the best-fit primordial power spectrum with a running from WMAP nine-year data with ns = 1.009 and αs = − 0.019. In panel b), we plot the LGMCA WMAP nine-year pseudo power spectrum (blue crosses) and the estimated instrumental noise power spectrum including the point sources power spectrum is shown (solid black line). The very small blue region corresponds to the 1σ interval around the mean reconstructed spectrum (i.e. blue region in panel a)). As can be seen, we do not detect a significant deviation of the WMAP nine-year data from the best fit near scale-invariant spectrum.

4. Conclusions

The primordial power spectrum describes the initial perturbations in the Universe and so provides an indirect probe of inflation or other structure-formation mechanisms. The simplest models of inflation are the most favoured by the data and predict a near scale-invariant power spectrum with a small running. One way to measure this spectrum is through the windows of the CMB data. The problem, however, is that the singular nature of the radiation transfer function and the joint estimation of the cosmological parameters in the transfer function and the primordial power spectrum, along with the different types of noise sources impose a limit into the full recovery of the primordial spectrum. Therefore, devising a technique that is sensitive enough to detect deviations from scale invariance is the key to recovering an accurate primordial power spectrum.

In this paper we have introduced a new non-parametric technique, called PRISM, to recover the primordial power spectrum from masked noisy CMB data. This is a sparse recovery method, which uses the sparsity of the primordial power spectrum as well as an adapted modelling for the noise of the CMB power spectrum. This algorithm assumes no prior shape for the primordial spectrum and does not require a coarse binning of the power spectrum, making it sensitive to both global smooth features (e.g. running of the spectral index) as well as local sharp features (e.g. a bump or an oscillatory feature). Another advantage of this method is that, thanks to the clever modelling of the sample variance on the input angular power spectrum, the regularisation parameter can be specified in terms of a signal-to-noise significance level for the detection of features. These advantages make this technique very suitable for investigating different types of departures from scale invariance in the primordial power spectrum, whether it is the running of the spectral index or some localised sharp features as predicted by some of the inflationary models.

We have investigated the strength of our proposed algorithm on a set of WMAP nine-year simulated data for three types of primordial power spectra: a near scale-invariant spectrum, a spectrum with a small running of the spectral index, and a spectrum with a localised feature. We have shown that our algorithm can easily recover the three spectra with an excellent accuracy in the range k ~ 0.001−0.1Mpc^-1. In addition, the errors in the recovered spectra are small enough that the three types of primordial spectra can easily be distinguished in the range k ~ 0.015−0.1Mpc^-1. This technique has proved that it can easily detect small global and localised deviations from a pure scale-invariant power spectrum and that it is suitable for distinguishing between simple models of the inflation.

Using PRISM, we have reconstructed a primordial power spectrum from the LGMCA WMAP nine-year data and have investigated possible departures from the WMAP nine-year near scale-invariant spectrum. We have not detected any significant deviations from this simple model of the primordial power spectrum. We have demonstrated the feasibility of using PRISM on masked CMB data contaminated by instrumental noise. Better constraints will be obtained in future works by processing Planck data which provides a much lower instrumental noise, thus improving the range of scales we are able to probe with much better accuracy.

To this end, we also acknowledge previous algorithms aimed at reconstructing the primordial power spectrum with no need for binning, most of which have been referenced in this paper. The most recent work is by Hunt & Sarkar (2014), who use the Tikhonov regularisation to solve the ill-posed inverse problem. The advantage of this approach is the linear relationship between the data and estimated primordial power spectrum. Just like our algorithm, the regularisation parameter has to be determined using simulations. However, because of the non-stationary noise in the data it is difficult to obtain a reliable regularisation parameter, whereas the robust noise handling in PRISM makes it easier to select an appropriate regularisation parameter. Furthermore, the wavelet regularisation in PRISM has the advantage of preserving the detected features well, which is not the case with Tikhonov regularisation. Another interesting work is by Hazra et al. (2013), who use an adapted and improved Richardson-Lucy algorithm, dubbed MRL, to reconstruct the primordial power spectrum. Because of the very high level of the instrumental noise on small scales in the WMAP nine-year data, the MRL algorithm takes the unbinned CMB spectrum only for ℓ< 900. For larger angular scales, ℓ = 900–1200, a binned CMB spectrum is used. In addition, because of the possible induced artefacts in the reconstructed primordial spectrum, a smoothing step may be necessary after the reconstruction is performed. Henceforth, compared to the MRL algorithm, the advantage of our algorithm is twofold. One is the ability to use the unbinned CMB spectrum for the whole multipole range ℓ = 2–1200 because of our accurate noise modelling on the CMB power spectrum. The other is that there is no need for an additional smoothing step in PRISM, as any residual features in the reconstructions are statistically significant.

The developed C++ and IDL codes will be released with the next version of iSAP (Interactive Sparse astronomical data Analysis Packages) via the website http://cosmostat.org/isap.html . All results have been obtained using the isap routine mrs_prism with the following command line:

pk = mrs_prism(Cl, noise=Nl, TransferMat=Mat)where Cl contains the observed pseudo-power spectrum of the masked noisy CMB maps, Nl is an estimate of the instrumental noise power spectrum N_ℓ, and Mat is the input transfer matrix which includes the effects of the radiation transfer function, the mask, the beam, and the HEALPix window (Eq. (39)). The transfer matrix can be computed using the isap routine mrs_transfer_matrix and by default the transfer matrix is derived from the WMAP nine-year best-fit cosmology model.

Acknowledgments

The authors would like to thank Amir Hajian, Gabriel Rilling, and Jeremy Rapin for their useful discussions. This work is supported by European Research Council grant SparseAstro (ERC-228261).

References

All Figures

	Fig. 1 A simulated noisy CMB map at 15 arcmin resolution obtained from LGMCA and masked with the WMAP kq85 mask. The noise level corresponds to the WMAP nine-year data. This map was generated from a CMB power spectrum for a primordial spectrum with ns = 0.972 and αs = 0.
In the text

Fig. 2 CMB pseudo power spectra for the three types of primordial power spectra. The solid blue line shows the pseudo spectrum based on a primordial spectrum with ns = 0.972 and αs = 0. The light blue crosses show one simulation of this spectrum, computed from the map in Fig. 1. The red line shows the pseudo spectrum for a primordial spectrum with ns = 0.972 and αs = − 0.017 and the orange line corresponds to a power spectrum with a localised feature at k = 0.03 Mpc-1. These spectra include the effects of the mask, the 15 arcmin beam, the HEALPix window for nside of 1024, and the instrumental noise power spectrum, which is shown by a solid black line.

In the text

Fig. 3 Reconstructions for the primordial power spectra and their corresponding CMB pseudo spectra. In blue we show the 2000 reconstructed spectra with ns = 0.972 and αs = 0 and in cyan the reconstruction for ns = 0.972 and αs = −0.017. In both cases the mean of the reconstructions is shown in orange and the fiducial input spectrum is shown in red. As can be seen, for k> 0.015Mpc-1 PRISM can reconstruct the primordial power spectra with such accuracy that the two are easily distinguishable, despite their very similar forms in Cℓ space. The shaded regions in panel b) correspond to the 1σ sample (cosmic) variance, which demonstrates the similarity of the two types of CMB spectra. The quality of the reconstruction can also be seen in the reconstructed angular power spectra which are extremely close to the theory and well within the 1σ sample variance intervals.

In the text

	Fig. 4 Reconstruction of the primordial power spectrum with ns = 0.972, αs = 0.0, and an additional feature around k = 0.03 Mpc-1 shown in green. The 2000 reconstructions are superimposed with their mean shown in orange. The fiducial input spectrum is shown in red. As can be seen, PRISM is able to recover both the position and the amplitude of the feature with great accuracy.
In the text

Fig. 5 Reconstruction of the primordial power spectrum from the LGMCA WMAP nine-year data and its corresponding pseudo spectrum shown in red. For comparison, we also show the mean of the reconstruction for ns = 0.972 and αs = 0 with a solid dark blue line with the 1σ interval around the mean shown as a shaded blue region. The WMAP nine-year fiducial primordial power spectrum with ns = 0.972 and αs = 0 is shown in yellow and in cyan we show the best-fit primordial power spectrum with a running from WMAP nine-year data with ns = 1.009 and αs = − 0.019. In panel b), we plot the LGMCA WMAP nine-year pseudo power spectrum (blue crosses) and the estimated instrumental noise power spectrum including the point sources power spectrum is shown (solid black line). The very small blue region corresponds to the 1σ interval around the mean reconstructed spectrum (i.e. blue region in panel a)). As can be seen, we do not detect a significant deviation of the WMAP nine-year data from the best fit near scale-invariant spectrum.

In the text