During the SO telecon of 2022-02-10 Jon and I mentioned that filter+bin could be expected to be just as vulnerable to loss of power at low l due to model errors as maximum-likelihood is. There was some skepticism about this claim, so here's a toy example demonstrating it. There are many ways model errors can cause bias, and in ACT we think we're dominated by unknown model errors. But for this simple demonstration we will consider sub-pixel errors, which arise from the true sky being smooth while we model it as pixellated. This will be enough to cause a loss of signal, but I want to emphasize that eliminating sub-pixel issues will just remove one of several possible sources of biases like this.
We consider a simple 1D case where 100 samples scan uniformly across 10 pixels:
npix = 10
nsamp = 100
pix = np.arange(nsamp).astype(float)*npix/nsamp
We use a standard nearest-neighbor pointing matrix, which means that each sample is assigned to the nearest pixel. Effectively each pixel is treated as a small rectangular bucket, and it doesn't matter where inside of a bucket a sample hits. This is the standard pointing matrix everybody uses when making sky maps.
P = np.zeros((nsamp,npix))
for i, p in enumerate(pix):
P[i,int(np.round(pix[i]))%npix] = 1
We assume a typical 1/f noise model. For the case of maximum-likelihood map-making, this will be our noise model, and for filter+bin it will be our filter. In either case we represent it as the matrix F.
freq = rfftfreq(nsamp)
nfreq= len(freq)
fknee= 0.03
# the max stuff just avoids making the mean infinitely uncertain
inv_ps = 1/(1+(np.maximum(freq,freq[1]/2)/fknee)**-3.5)
def filter_func(x):
return np.fft.irfft(inv_ps*np.fft.rfft(x), n=len(x))
F = func_to_mat(filter_func, nsamp)
Here's what the noise power spectrum looks like. The noise weighting/filter is just the inverse of this.
For the signal we choose a simple, large-scale sine wave. We don't add any noise to this because all the mapmaking algorithms are linear in the data, so adding noise would just make our maps more noisy without changing their mean value.
def build_signal(pix): return np.sin(2*np.pi*pix/npix)
signal = build_signal(pix)
Ok, are now ready to make some maps. The first map will be a simple, unweighted binned map.
map_binned = np.linalg.solve((P.T.dot(P)), P.T.dot(signal))
Next is a maximum-likelihood map.
map_ml = np.linalg.solve((P.T.dot(F).dot(P)),P.T.dot(F.dot(signal)))
We also make a filter+bin map (map_fb), which we then
debias using an observation matrix. This turns out to
be equivalent to maximum-likelihood mapmaking.
map_fb = np.linalg.solve(P.T.dot(P), P.T.dot(F).dot(signal))
obsmat = np.linalg.inv(P.T.dot(P)).dot(P.T.dot(F).dot(P))
map_fb_deobs = np.linalg.solve(obsmat, map_fb)
We can also use a suite of simulations to measure and deconvolve the filter+bin transfer function.
nsim = 1000
sim_ips = np.zeros(npix//2+1)
sim_ops = np.zeros(npix//2+1)
for i in range(nsim):
sim_imap = np.random.standard_normal(npix)
sim_omap = np.linalg.solve(P.T.dot(P), P.T.dot(F).dot(P).dot(sim_imap))
sim_ips += np.abs(np.fft.rfft(sim_imap))**2
sim_ops += np.abs(np.fft.rfft(sim_omap))**2
tf = (sim_ops/sim_ips)**0.5
map_fb_detrans = np.fft.irfft(np.fft.rfft(map_fb)/tf, n=npix)
Finally, we make a generalized destriped map in the limit of one baseline per sample, which should be equivalent to ML.
I = np.eye(nsamp)
iCa = np.linalg.inv(np.linalg.inv(F) - I)
Z = I-P.dot(np.linalg.solve(P.T.dot(P), P.T))
a = np.linalg.solve(Z+iCa, Z.dot(signal))
map_ds = np.linalg.solve(P.T.dot(P), P.T.dot(signal - a))
Here's what these 1D "maps" look like after deconvolving the pixel window, compared to the signal itself.
As you can see, all methods except for plain binning result in a strongly biased map, and filter+bin is just as biased as maximum likelihood. Another interesting reslut is that using simulations to debias filter+bin resulted in a worse residual bias than when debiasing using an obseration matrix. This might be due to mode mixing, which the transfer function ignores, unlike the observation matrix. If I add a big sine mode to the simulations, then the trasnfer function approach gets much better.
Finally, let's check what happens when the data actually does follow the model, just to see that my implementation of the equations isn't buggy.
signal2 = P.dot(map_binned)
map_binned2 = np.linalg.solve(P.T.dot(P), P.T.dot(signal2))
map_ml2 = np.linalg.solve((P.T.dot(F).dot(P)),P.T.dot(F.dot(signal2)))
map_fb2 = np.linalg.solve(P.T.dot(P), P.T.dot(F).dot(signal2))
map_fb_deobs2 = np.linalg.solve(obsmat, map_fb2)
map_fb_detrans2 = np.fft.irfft(np.fft.rfft(map_fb2)/tf, n=npix)
a2 = np.linalg.solve(Z+iCa, Z.dot(signal2))
map_ds2 = np.linalg.solve(P.T.dot(P), P.T.dot(signal2 - a2))
This results in the following graph:
Here all mapmaking methods exactly reproduce the signal in the pixel centers, except for filter+bin with transfun debiasing.
Anyway, I hope this demonstrates that model errors is something we should be very wary of in SO, regardless of what type of mapmaking we use! And again, even though I used subpixel errors to demonstrate this, there are many other types of model error that could have this type of effect, so just solving the subpixel issue won't solve the whole problem.
You can download the full program here, and an interactive jupyter notebook is here.