The reason we can is because the Planck satellite moves very slowly, so it measures the same point on the sky several times as it circles around. So, we have a set of data that has several measurements for each pixel for some number of pixels. In order to remove the drift we make the logical assumption that the true temperature is the one that minimizes the variance of the measurements over the entire sky, which we find by minimizing the variance in each individual pixel. It does this by finding the coeffecients of smooth step functions, which are generated by splines. A note about the matrix inversion.
Alpha is, by definition, symmetric.
And since it is almost always positive definite, the quickest and most efficient algorithm is Cholesky Decomposition. When this is the case, the algorithm uses the slower but much more stable matrix inversion routine using Singular Value Decomposition. Here is a short description of each of the subroutines and their purpose in the program. This function creates our original map and simulates the kind of data we will be getting from Planck. Currently, it simulates a map of a sky with a constant underlying value, with Gaussian noise and linear and sinusoidal drift added to that.