Cooling Improves Cosmic Microwave Background Map-Making When Low-Frequency Noise is Large

In the context of Cosmic Microwave Background data analysis, we study the solution to the equation that transforms scanning data into a map. As originally suggested in"messenger"methods for solving linear systems, we split the noise covariance into uniform and non-uniform parts and adjust their relative weights during the iterative solution. With simulations, we study mock instrumental data with different noise properties, and find that this"cooling"or perturbative approach is particularly effective when there is significant low-frequency noise in the timestream. In such cases, a conjugate gradient algorithm applied to this modified system converges faster and to a higher fidelity solution than the standard conjugate gradient approach. We give an analytic estimate for the parameter that controls how gradually the linear system should change during the course of the solution.


INTRODUCTION
In observations of the Cosmic Microwave Background (CMB), map-making is an intermediate step between the collection of raw scanning data and the scientific analyses, such as the estimation of power spectra and cosmological parameters. Next generation CMB observations will generate much more data than those today, and so it is worth exploring efficient ways to process the data even though, on paper, the map-making problem has long been solved.
The time-ordered scanning data is summarized by where d, m, and n are the vectors of time-ordered data (TOD), the CMB sky-map signal, and measurement noise. P is a sparse matrix in the time-by-pixel domain that encodes the telescope's pointing. Of several map-making methods (Tegmark 1997), one of the most common is the method introduced for the Cosmic Background Explorer (COBE, Janssen & Gulkis 1992). This optimal, linear solution is wherem provides the standard generalized least squares minimization of the χ 2 statistic, Here we assume that the noise has zero mean n = 0, and the noise covariance matrix N = nn † is diagonal in frequency space. In the case where the noise is Gaussian, the COBE solution is also the maximum likelihood solution.
With current computational power, we cannot solve form by calculating P † N −1 P −1 P † N −1 d directly.
The noise covariance matrix N is often sparse in the frequency domain and the pointing matrix P is sparse in the time-by-pixel domain. In experiments currently under design, there may be ∼ 10 16 time samples and ∼ 10 9 pixels, so these matrix inversions are intractable unless the covariance is uniform (proportional to the identity matrix I). We can use iterative methods, such as conjugate gradient descent, to avoid the matrix inversions, and execute each matrix multiplication in a basis where the matrix is sparse, using a fast Fourier transform to go between the frequency and time domain.
As an alternative to conjugate gradient descent, Huffenberger & Naess (2018) showed that the "messenger" iterative method could be adapted to solve the linear map-making system, based on the approach from Elsner & Wandelt (2013) to solve the linear Wiener filter. This technique splits the noise covariance into a uniform part and the remainder, and introduces an additional vector that represent the signal plus uniform noise. This messenger field acts as an intermediary between the signal and the data and has a covariance that is conveniently sparse in every basis. Elsner & Wandelt (2013) also introduced a cooling scheme that takes advantage of the split covariance: over the course of the iterative solution, we adjust the relative weight of those two parts. Starting with the uniform covariance, the modified linear system gradually transforms to the final system, under the control of a cooling parameter. In numerical experiments, Huffenberger & Naess (2018) found that a map produced by the cooled messenger method converged significantly faster than for standard conjugate gradient methods, and to higher fidelity, especially on large scales. Papež et al. (2018) showed that the the messenger field approach is equivalent to a fixed point iteration scheme, and studied its convergence properties in detail. Furthermore, they showed that the split covariance and the modified system that incorporates the cooling can be solved by other means, including a conjugate gradient technique, which should generally show better convergence properties than the fixed-point scheme.
With simulations, we have confirmed this conclusion. In their numerical tests, Papež et al. (2018) did not find benefits to the cooling modification of the map-making system, in contrast to the findings of Huffenberger & Naess (2018).
In this paper, we show that the difference arose because the numerical tests in Papež et al. (2018) used much less low-frequency (or 1/f ) noise than Huffenberger & Naess (2018), and show that the cooling technique improves map-making performance especially when the low-frequency noise is large. This performance boost depends on a proper choice for the pace of cooling. Kodi Ramanah et al. (2017) showed that for Wiener filter the cooling parameter should be chosen as a geometric series. In this work, we give an alternative interpretation of the parameterizing process and show that for map-making the optimal choice (unsurprisingly) is also a geometric series.
In Section 2 we describe our methods for treating the map-making equation and our numerical experiments. In Section 3 we present our results. In Section 4, we list our conclusions. In Appendix A we derive the prescription for our cooling schedule.

Parameterized Conjugate Gradient Method
The messenger field approach introduced an extra cooling parameter λ to the map-making equation, and solved the linear system with an alternative parameterized covariance N (λ) = λτ I +N . The parameter τ = min(diag(N )) represents the uniform level of (white) noise in the original covariance. The remain-derN ≡ N − τ I is the non-uniform part of the original noise covariance. (Here N without any arguments denotes the original noise covariance matrix N = nn † .) In this work we find it more convenient to work with the reciprocal of cooling parameter η = λ −1 which represents the degree of heteroscedasticity (non-uniformity) in the parameterized covariance When η = 1 this parameterized covariance N (η) equals N . Papež et al. (2018) showed that the conjugate gradient method can be easily applied to the cooled map-making problem. In our notation, this is equivalent to iterating on the parameterized map-making equation as we adjust the parameter through a set of levels {η i }.
(We usem with no η argument to mean the estimatedm in Eq. 2, independent of η.) This equation leads to the same system as the paramaterized equation in the messenger field method, because N (η) = λ −1 N (λ) and the condition number do not change upon scalar multiplication to both sides of the equation. For concreteness we fix the preconditioner to M = P † P for all calculations. When η = 0, the noise covariance matrix N (0) is homoscedastic (uniform), and the solution is given by the simple binned mapm(0) = P † P −1 P † d, which can be solved directly.
Since the non-white partN is the troublesome portion of the covariance, we can think of the η parameter as increasing the heteroscedasticity of the system, adding a perturbation to the solution achieved at a particular stage, building ultimately upon the initial uniform covariance model. Therefore, this quasi-static process requires η increase as 0 = η 0 ≤ η 1 ≤ · · · ≤ η final = 1, at which point we arrive at the desired map-making equation, and the solutionm(1) =m.
We may iterate more than once at each intermediate η i : we solve equation (5) with conjugate gradient iterations using the result from the previous calculation m(η i−1 ) as the initial value. We move to next parameter η i+1 when the norm of residual vector is an order of magnitude smaller than the norm of the right hand side of Eq. 5.
This is not stringent enough to completely converge at this η-level, but we find that it causes the system to converge sufficiently to allow us to move on to the next η.

Analytical expression for {η i } series
The next question is how to appropriately choose these monotonically increasing parameters η. We also want to determine η 1 , · · · , η n−1 before starting conjugate gradient iterations, because the time ordered data d is very large, and we do not want to keep it in the system memory during calculation or repeatedly read it in from disk. If we determine η 1 , · · · , η n−1 before the iterations, then we can precompute the right-hand side of Eq. 5 for each η i and keep these map-sized objects, instead of the entire time-ordered data.
In Appendix A, we show that a generic good choice for the η parameters is given by this geometric series whereN f are the eigenvalues ofN under frequency representation. This is one of our main results. It not only tells us how to choose parameters η i , but also when we should stop the perturbation, and set η = 1. For example, if the noise covariance matrix N is almost uniform, thenN = N − τ I ≈ 0, and we would have τ /max(N f ) > 1. This tell us that we don't need to use the parameterized method at all, because η 0 = 0 and η 1 = η 2 = · · · = 1. This corresponds to the standard conjugate gradient method with simple binned map as the initial guess (as recommended by Papež et al. 2018).

Intuitive Interpretation of η
Here is a way to interpret the role of η that is less technical than Appendix A. Our ultimate goal is to find m(1) which minimizes χ 2 (m) in Eq. 3. Since N is diagonal in frequency space, χ 2 could be written as a sum of all frequency modes | The weight is large for low-noise frequency modes (small N f ), and small for high-noise modes. Which means χ 2 (m) would favor the low-noise modes, and therefore the conjugate gradient map-making focuses on minimizing the error ε ≡ d − P m in the low-noise part.
After introducing η, we minimize χ 2 (m, η) in Eq. A1 instead. For η = 0, N −1 (0) ∝ I the system is homoscedastic and the estimated mapm(0) does not prioritize any frequency modes. As we slowly increase η, we decrease the weight for the high-noise modes, and focusing minimizing error for the low-noise part. If we start with η 1 = 1 directly, which corresponds to the vanilla conjugate gradient method, then the algorithm will focus most on minimizing the low-noise part, such that χ 2 would converge very fast on the low-noise modes (typically high temporal frequencies and small spatial scales), but slowly on the high-noise modes (low frequencies and large scales). However by introducing the η parameter, we let the solver first treat every frequency equally. Then as η slowly increases, it gradually gives more focus to the lowest noise part.

Computational Cost
To properly compare the performance cost of this method with respect to the vanilla conjugate gradient method with the simple preconditioner, we need to compare their computational cost at each iteration.
We could define A(η) ≡ P † N (η) −1 P and b(η) ≡ P † N (η) −1 d, and equation 5 could be written as A(η i )m(η i ) = b(η i ). The right-hand side b(η i ) could be computed before iterating, since we have determined {η i } in advance, so it will not introduce extra computational cost. The most demanding part of conjugate gradient method is calculating its left hand side A(η i )m, because it contains a Fourier transform of P m from the time domain to frequency domain and an inverse Fourier transform of N (η i ) −1 P m from the frequency domain back to time domain, which is order O(n log n) with n being the length of time ordered data. Compared to the traditional conjugate gradient method, we swap N −1 with N (η) −1 , and the cost is the same for one step, since both methods need a fast Fourier transform and inverse fast Fourier transform at one iteration.
At each η i level, we use the residual to determine whether to switch to the next level (η i+1 ), as is Equation (7). Calculation of the residual vector r(m, η i ) is part of the conjugate gradient algorithm, so this will not add extra cost either. Therefore, overall introducing the η will not have extra computational cost within the conjugate gradient iterations.
However, we start a new conjugate gradient algorithm whenever η i updates to η i+1 . Thus we must re-initialize the conjugate gradient algorithm, re-calculating the residual r(m, η i+1 ) based on new η i+1 . This residual calculation contains an extra A(η i )m operation. Therefore, if we have a series η 1 , η 2 , η 3 , · · · , η nη , there will have n η − 1 extra A(η)m operations compare to the traditional conjugate gradient method. If the total number of iterations is much larger than n η , then this extra cost is negligible. For our simulation, this extra step would have rather significant impact on final result. To have a fair comparison between the parameterized and traditional conjugate gradient method, we will present our results with number of P † N (η) −1 P m operations as horizontal axis.

Numerical Simulations
To compare these algorithms, we need to do some simple simulations of scanning processes, and generate the Noise power spectra that we use in our mapmaking simulations. These show a variety of low-frequency behavior, parameterized by Eq. 9, with white noise at high frequency and a low-frequency power-law slope α = 3. Here we show two knee frequencies, f knee = 10 Hz (solid lines) and f knee = 0.1 Hz (dashed lines). For each knee frequency, we have shown an unflattened spectrum (fapo = 0 Hz), and two flattened ones (fapo = 0.1f knee and 0.01f knee ). The vertical line shows our scanning frequency.
time ordered data from a random sky signal. 1 Our sky is a small rectangular area, with two orthogonal directions x and y, both with range from −1 • to +1 • . The signal has Stokes parameters (I, Q, U ) for intensity and linear polarization. For the scanning process, our mock telescope contains nine detectors, each with different sensitivity to polarization Q and U . It scans the sky with a raster scanning pattern. Its scanning frequency is f scan = 0.1 Hz and sampling frequency is f sample = 100 Hz. The telescope scans the sky horizontally then vertically. This gives the noiseless signal s. The sky signal in the timestream has a root-mean-square (RMS) of 56 µK. The signal is continuous, so that it has structure on sub-pixel scales, but we find that our main conclusions remain the same when the input signal is pixelized. In map-making, we digitize the position (x, y) into 512 × 512 pixels.
We model the noise power spectrum with The source code and other information are available at https: //github.com/Bai-Chiang/CMB map making with cooling which is white at high frequencies, a power law below the knee frequency, and gives us the option to flatten the low-frequency noise below an apodization frequency (like in Papež et al. 2018). Note that as f apo → 0, P (f ) → σ 2 (1 + (f /f knee ) −α ), and it becomes a 1/f -type noise model. Dünner et al. (2013) measured the slopes of the atmospheric noise in the Atacama under different water vapor conditions, finding α = 2.7 to 2.9. Here we use σ 2 = 10 µK 2 , α = 3, and compare the performance under different noise models. In our calculations, we choose different combinations of f knee and f apo as in Figure 1. The noise spectra with the most low frequency noise have high f knee or low cut-off f apo .
The noise covariance matrix is a diagonal matrix in frequency space, where ∆ f is equal to the reciprocal of total scanning time T ≈ 1.05× 10 4 seconds. Finally, we get the simulated time ordered data d = s + n by adding up the signal and noise.

RESULTS
We compare the standard conjugate gradient method versus the conjugate gradient with our perturbed linear system. Both methods use the simple preconditioner P † P . Figure 2 shows the χ 2 results for 1/f noise models (f apo = 0) with different knee frequencies. Note that the χ 2 values in all figures are calculated based on the standard χ 2 (m) in Eq. 3, not the η-dependent χ 2 (m, η) of the modified system (Eq. A1). The minimum χ 2 min that we use for comparison is calculated from a deliberately slowed and well-converged parameterized conjugate gradient method: one with 100 η values and that halts when the final norm of the residual r(m, 1) is smaller than 10 −5 × P † N −1 d , or 100 iterations after η = 1. From Figure 2, we can see for the 1/f noise model, when f knee 10f scan the parameterized method starts showing an advantage over the vanilla conjugate gradient method.
In Figure 3, we fixed f knee = 10 Hz, and change f apo . As we decrease f apo relative to f knee , increasing the amount of low-frequency noise, the parameterized conjugate gradient method performs better.
Looking at the power spectrum in Figure 1, when f knee is small, or f apo is large, there is not much low-frequency noise. These situations corrrespond to the left-side plots in Figure 2 and Figure 3. The right-side graphs have significant amount of low frequency noise. We conclude that the introduction of the slowly-varying η parameter (f knee , f apo ) = (1, 0) Hz CG CG with Figure 2. Convergence properties depend on the amount of low-frequency noise, which increases from the left panel to the right panel with increasing knee frequency. The map-making equation 2 minimize the χ 2 (m), so the curve which falls fastest versus the number of operations is the preferred method. We compare the traditional conjugate gradient method ("CG," blue line) with the parameterized conjugate gradient method ("CG with η," orange line) under different 1/f noise models (fixed fapo = 0 Hz but different f knee in Eq. 9). When f knee 10 fscan = 1 Hz, the significant amount of low-frequency noise causes the parameterized conjugate gradient method to start showing its advantage. The vertical axis is rescaled such that all curves start from 1.   Figure 2, low-frequency noise increases from left to right, but by flattening the low-frequency noise at an apodization frequency. Low-frequency noise increases with decreasing apodization frequency (compare Figure 1). We again compare the traditional conjugate gradient method ("CG," blue line) with the parameterized conjugate gradient method ("CG with η," orange line). When fapo is much smaller than f knee , there is a lot of low-frequency noise and the parameterized conjugate gradient method is better (ultimately falls faster) than the traditional one.
improves performance most when there are large lowfrequency noise contributions. We also tried different 1/f noise slopes α. For α = 2, the conclusion is the same as α = 3. When α = 1, the low-frequency noise is reduced compared to the cases with steeper slopes, and the vanilla conjugate gradient method is preferred, except some cases with very large knee frequency like f knee = 100 Hz and f apo = 0 which favors the parameterized method. In Papež et al. (2018), the slope α = 1 and the noise power spectrum is flattened at f apo ≈ 0.1f knee . Their knee frequency is the same as their scanning frequency, so is most like our case when f knee = f scan = 0.1 Hz. Their case had little low-frequency noise, and we confirm their specific result that the standard conjugate gradient method converges faster in that case. In general, however, we find cases with significantly more low-frequency noise benefit from the cooling/parameterized approach.

CONCLUSIONS
We analyzed the parameterized conjugate gradient map-making method that is inspired by the messengerfield idea of separating the white noise out of the noise covariance matrix. Then we gave an analytical expression for the series of η parameters that govern how quickly the modified covariance adjusts to the correct covariance, and showed that this method adds only the extra computational cost of re-initializing the conjugate gradient process based on a new η parameter.
We tested this method for different noise power spectra, both flattened and non-flattened at low frequency. The results showed that the parameterized method is faster than the traditional conjugate gradient method when there is a significant amount of low-frequency noise. It could be further improved if we could get a more accurate estimation for the change in χ 2 as a function of the η parameter, either before iteration or without using time ordered data during iteration.
Also note that we fixed the preconditioner as M = P † P during our calculation, this parameterizing process could be applied to any preconditioner and possibly improve performance when there is significant amount of low-frequency noise.
This type of analysis for the cooling parameter may also be apply to other areas, like the Wiener filter. Papež et al. (2018) showed that the messenger field method of Elsner & Wandelt (2013) for solving Wiener filter problem could also be written as a parameterized conjugate gradient algorithm. It stands to reason that such a system may also benefit from splitting and parameterizing its noise covariance, depending on the noise properties. The benefits to map-making from a cooled messenger method seem to come from the cooling and not actually from the messenger field that inspired it. However, the messenger field approach may still have a role in the production and analysis of CMB maps. In particular, the close connection between the messenger method and Gibbs sampling may allow us to cheaply generate noise realizations of a converged map by generating samples from the map posterior distribution, something that we will continue to explore in future work.