-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathindex.html
More file actions
1270 lines (1232 loc) · 62.5 KB
/
index.html
File metadata and controls
1270 lines (1232 loc) · 62.5 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
<html>
<head>
<title>How to avoid X'es around point sources in maximum likelihood CMB maps</title>
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
extensions: ["tex2jax.js"],
jax: ["input/TeX", "output/HTML-CSS"],
tex2jax: {inlineMath: [["$","$"],["\\(","\\)"]]}
});
</script>
<script type="text/javascript" src="http://folk.uio.no/sigurdkn/mathjax/MathJax.js"></script>
<style>
body { max-width: 1200px; margin-left: auto; margin-right: auto; padding-left: 1em; padding-right: 1em; }
img { width: 100%; }
figure { display: block; margin-left: auto; margin-right: auto; width: 100%; max-width: 40em; border: dashed black 1px; margin: 0.5em; }
figure.wide { max-width: 75em; }
table { table-layout: fixed; width: 100%; }
table.pix img { image-rendering: optimizeSpeed; }
table.text td { text-align: center; }
table.fit { table-layout: auto; }
td.good { font-weight: bold; color: green; }
td.medi { font-weight: bold; color: orange; }
td.bad { font-weight: bold; color: red; }
td.neutral { }
td.mixed { font-weight: bold; color: #05a; }
h1, h2, h3, h4, h5 { text-align: center; }
p, figcaption, li, .abstract { text-align: justify; }
figcaption { margin: 0.5em; }
tr.sub th { font-size: 70%; }
dfn { font-weight: bold; }
.authlist, .affil li { text-align: center; }
.author { font-size: 120%; }
.afftag { vertical-align: super; font-size: 70%; }
.affil { list-style: none; }
div.abstract { margin-left: 5em; margin-right: 5em; font-size: 90%; }
</style>
</head>
<body>
<h1>How to avoid X'es around point sources in maximum likelihood CMB maps</h1>
<div class=authlist>
<span class=author>Sigurd K. Næss</span><span class=afftag>1</span>
</div>
<ul class=affil>
<li><span class=afftag>1</span> Center for Computational Astrophysics, Flatiron Institute</li>
</ul>
<div class=abstract>
The maximum likelihood estimator for CMB map-making is optimal and unbiased as long as
the data model is correct, but in practice it rarely is, with model errors including
sub-pixel structure and instrumental problems like time-variable gain and pointing errors.
In the presence of such errors, the solution is biased, with the local error in each pixel
leaking outwards along the scanning pattern by a noise correlation length. The most important
sources of such leakage are strong point sources, and for common scanning patterns the leakage
manifests as an X around each such source. I discuss why this happens, and present several
old and new methods for mitigating and/or eliminating this leakage, along with a small
stand-alone TOD simulator and map-maker in Python that implements them.
</div>
<h2>1. Introduction</h2>
<p>Cosmic Microwave Background (CMB) telescopes map the sky by scanning an array of detectors repeatedly across
the sky, resulting in a time-series of samples $d$. This process is usually
modeled via the linear equation (<a href=https://arxiv.org/abs/astro-ph/9611130>Tegmark 1997</a>)</p>
\begin{align}
d &= Pm + n \label{eq:model} \tag{1}
\end{align}
<p>where $d$ is the vector of samples,
$m$ is a vector representing the pixels of the sky we are trying to reconstruct, and
$P$ is a response matrix that encodes where the telescope was pointing on the sky
when each sample was taken and how it responded to that. $n$ represents the
noise contribution to each sample, and is usually taken to be normal distributed
with some covariance $N$. The maximum likelihood solution to this equation is </p>
\begin{align}
\hat m &= (P^T N^{-1} P)^{-1} P^T N^{-1} d \label{eq:mapmaking} \tag{2}
\end{align}
<p>and it's simple to see that this solution is unbiased by inserting our model for $d$:</p>
\begin{align}
\langle \hat m \rangle &= (P^T N^{-1} P)^{-1}P^T N^{-1}(Pm + \overbrace{\langle n \rangle}^{0}) = m
\end{align}
<p>Given this result it might be surprising to hear that most maximum-likelihood maps
made using eq. \ref{eq:mapmaking} are biased to some extent, with the archetypical example
being lines extending away from bright point sources in the map, often in the shape of an X
(but this depends on the telescope's scanning pattern). The source of this bias is <em>model
error</em>: the failure of eq. \ref{eq:model} to accurately describe the real data; and
in this paper I will describe the main types of model error in maximum-likelihood mapmaking, simulate their
effect and investigate several methods for mitigating or eliminating them.</p>
<h2>2. Sub-pixel structure</h2>
<p>The most obvious problem with eq. \ref{eq:model} is that it models the sky as a
finite vector of pixels, which will leave out any signal on scales smaller than the
pixel spacing. This is exacerbated by the canonical choice of a <em>nearest-neighbor</em> response
model in $P$.
Every CMB map-maker I am aware of uses nearest-neighbor, including destripers like
MADAM (<a href=https://arxiv.org/abs/0907.0367>Keihänen et al. 2010</a>),
and SRoll (<a href=https://arxiv.org/abs/1807.06207>Planck Collaboration, 2018</a>),
filter+bin map-makers like those used in SPT (<a href=https://arxiv.org/abs/1111.7245>Schaffer et al. 2011</a>)
and BICEP (<a href=https://arxiv.org/abs/1403.3985>BICEP2 Collaboration, 2014</a>),
or the maximum likelihood map-makers
used in ACT (<a href=https://arxiv.org/abs/1610.02360>ACT Collaboration, 2017</a>)
and QUIET (<a href=https://arxiv.org/abs/1508.02778>Ruud et al. 2015</a>).
In this model, the value of each sample is simply given by the value of the
closest pixel to it, regardless of how far away from the pixel center it is. This results
in a response matrix with a very simple structure: For each row (representing a sample)
a single element will be 1 (representing the pixel hit by that sample), and all others
will be zero.</p>
<p>Multiplication by a matrix with this structure is simple and efficient.
Given a map m[ny,nx] and arrays y[nsamp], x[nsamp] containing the x and y pixel coordinates
of each sample, the forward operation $d = Pm$ can be implemented as</p>
<pre><code>for i in range(nsamp):
d[i] = m[round(y[i]),round(x[i])]</code></pre>
<p>and the transpose operation $m = P^Td$ can be implemented as</p>
<pre><code>for i in range(nsamp):
m[round(y[i]),round(x[i])] += d[i]</code></pre>
<p>However, this simplicity comes at a cost: The signal model implied by this
response matrix has a constant value inside each pixel and a discontinuous jump
at the edge of each pixel. This model and the resulting error is illustrated in figure 1.
</p>
<figure class=wide>
<table>
<tr><th>uncorrelated</th><th>correlated</th></tr>
<tr>
<td><img src=examples/src_1d_uncorr_tod_log.svg></td>
<td><img src=examples/src_1d_corr_tod_log.svg></td>
</tr>
</table>
<figcaption>
<dfn>Figure 1</dfn>: A simulated data vector consisting of a noiseless scan through a smooth, Gaussian profile,
compared to its nearest-neighbor-pixellated maximum likelihood
model. <dfn>Left</dfn>: When an uncorrelated noise model is used ($N \propto I$)
the residual is large but stays in the pixel where they were sourced. This error can
be characterized as a simple pixel window. <dfn>Right</dfn>: When a correlated noise
model is used the residuals are still large, but are now spread out over a larger area
than the signal that sourced them.
</figcaption>
</figure>
<p>A step-function model like this clearly cannot match the smooth behavior of the beam-convolved
sky, and the only term in eq. \ref{eq:model} that can absorb such a mismatch between the data
and the model is the noise term $n$. How this model error manifests is therefore determined by
our noise model $N$.</p>
<p>If $N$ is assumed to be uncorrelated ($N$ diagonal in time domain), then a high value of the residual
(honorary noise as far as our data model is concerned) in one sample does not tell us
anything about what the noise is doing in any other samples, and so the residual
just stays in the pixel that birthed it. The value in each pixel is therefore just
the mean of the samples that hit it (left panel of figure 1).
This is approximately
the mean value of the signal inside the area covered by the pixel, which is a straightforward
and useful thing to measure, and is not the sort of bias we are worried about here.</p>
<p>If, on the other hand, $N$ is assumed to be correlated, then the presence of high
"noise" (actually residual) in some samples
will make the maximum likelihood prection for other samples within a noise correlation length
have a similar value for the noise. For long correlation lengths this could affect samples
many pixels away from the source of the residual. This is illustrated in the right panel
of figure 1. Here large sub-pixel residuals in the center result in a non-zero expectation
value for the noise at large distances. Given that the actual data has negligibly small values
there, the best fit model must cancel the expected noise, resulting in a non-zero model
extending far beyond the area with significant signal. Effectively the correlated
noise model is causing local model errors to leak roughly a correlation length
into the surrounding pixels.</p>
<h3>2.1. Higher-order mapmaking</h3>
<p>Since the largest source of sub-pixel bias is the use of nearest-neighbor interpolation in
the response matrix $P$, an obvious next step is to use a smoother interpolation function,
such as bilinear and bicubic spline interpolation. These are popular in image processing, but
are as far as I am aware unused in CMB map making.</p>
<p>Bilinear interpolation is considerably slower and more complicated than nearest neighbor,
and each sample now interacts with the four closest pixels instead of just one.
The forward operation $d=Pm$ is now (based on the implementaiton of <code>scipy.ndimage.map_coordinates</code>):</p>
<pre><code>for i in range(nsamp):
py1, px1 = floor(y[i]), floor(x[i])
py2, px1 = py1+2, px1+1
ry, rx = y[i]-py1, x[i]-px1
vx1 = m[py1,px1] * (1-rx) + m[py1,px2] * rx
vx2 = m[py2,px1] * (1-rx) + m[py2,px2] * rx
d[i] = vx1 * (1-ry) + vx2 * ry</code></pre>
<p>From this we can derive the The transpose operation $m=P^Td$ is:</p>
<pre><code>for i in range(nsamp):
py1, px1 = floor(y[i]), floor(x[i])
py2, px1 = py1+2, px1+1
ry, rx = y[i]-py1, x[i]-px1
vx1 = d[i] * (1-ry)
vx2 = d[i] * ry
m[py1,px1] += vx1 * (1-rx)
m[py1,px2] += vx1 * rx
m[py2,px1] += vx2 * (1-rx)
m[py2,px2] += vx2 * rx</code></pre>
<p>Bicubic spline interpolation is another large step up in complexity. See the appendix for details.</p>
<p>Figure 2 shows how these higher-order interpolation schemes perform for the same
1D Gaussian simulation we considered in figure 1.</p>
<figure class=wide>
<img src=examples/src_1d_highorder.svg></td>
<figcaption>
<dfn>Figure 2</dfn>: Comparison of the maximum-likelihood models for nearest neighbor,
linear and cubic interpolation in the response matrix for uncorrelated and correlated
noise models, for the same data as in figure 1. Linear and cubic interpolation fit the
central part of the profile much better than nearest neighbor, leading to lower discrepancy
between the data and model at high radius in the correlated case (solid curves).
But this situation is reversed for the uncorrelated
noise model (dashed curves). This happens because the linear and cubic models are inherently non-local,
and are willing to sacrifice accuracy in areas where the signal is very weak in order to
get a better match in areas with strong signal.
</figcaption>
</figure>
<p>Both linear and cubic interpolation fit the peak of the Gaussian much better than nearest neighbor,
resulting in roughly 4 orders of magnitude lower spurious signal at high radius for the case with
a correlated noise model. However, they only reach this level at large radii. Before this, the
spurious signal is dominated by a new effect that was not present for nearest neighbor interpolation:
Higher order interpolation is inherently non-local, which means that each pixel is correlated with
its neighbors (and neighbors-neighbors for bicubic). This is responsible for the exponentially
decaying ringing pattern we see in the figure for these interpolation types, even for uncorrelated noise.
Effectively the fit is sacrificing accuracy at large radius in order to fit the strong, central
peak slightly better.</p>
<p>So we see that high interpolation order has smaller total residuals and hence lower
non-local artifacts from the noise correlations, but low-order interpolation has shorter
correlation length for the interpolation functions themselves, leading to a faster decay of
the ringing. For the example in figure 2, the sweet spot is linear interpolation.</p>
<p>Even better locality (at even higher cost) should be possible with an interpolation
function based on the actual point spread function of the instrument (cite Wandelt here?).
However, there are diminishing returns when reducing sub-pixel errors because there is another
class of model errors that if often much larger, as we shall see in section 3.</p>
<h4>2.1.1 Iterative approximation</h4>
<p>Since the response matrix often is the bottleneck in maximum-likelihood mapmaking,
it may be too expensive to make this step several times slower by using higher-order
mapmaking. A simple approximation that avoids most of this slowdown is to first make
a normal map using the fast, nearest-neighbor response matrix, and then subtract this
estimate of the sky from the time-ordered data using a bilinear or bicubic spline
response matrix. This process can be repeated if necessary, though in my tests there
was little point in doing more than 1-3 iterations.
\begin{align}
\hat m_0 &= 0 \\
\hat m_{i+1} &= \hat m_i + (P^TN^{-1}P)^{-1}P^TN^{-1}(d - \bar P \hat m_i) \label{eq:iterative} \tag{3}
\end{align}
Here $P$ is the standard nearest-neighbor response matrix and $\bar P$ is
the higher-order one we are approximating. Note that this iteration converges
on an <em>approximation</em> to the higher-order map, not the exact solution
we would have gotten by using $\bar P$ in the equation \ref{eq:mapmaking}.
However, as we shall see in the results section it does a pretty good job.</p>
<p>If one is already doing iterative mapmaking to avoid bias from
estimating the noise model $N$ from the full data $d$ rather than
just the noise $n$, then this iterative approximation to higher-order mapmaking
can simply be done as a part of that, and will incur almost no extra cost.</p>
<p>Alternatively, if one is solving for each map using an iterative method
like preconditioned conjugate gradients (CG), then the small scales that are the
main source of sub-pixel error often converge much faster than the larger ones.
In this case only a few CG steps are needed for all but the last iteration
of equation \ref{eq:iterative}, again resulting in almost no overhead compared to
standard mapmaking.</p>
<h3>2.2. Point source subtraction</h3>
<p>In figure 1 we saw that non-local sub-pixel artifacts could leak far away from the
signal that sourced them, but it also bears noting that their amplitude was about
$10^4$ times lower than the source. For a smooth, low-contrast signal like the
cosmic microwave background, each pixel is contaminated by leakage from
other pixels in the vicinity, but since those other pixels have signal of similar
amplitude but random phase, one can expect each pixel to only be contaminated at the
$10^{-4}$ level, which is negligible. Non-local model errors are therefore only
a concern when a strong, compact signal is close to a region of weaker signal. In
CMB map-making the typical case would be a strong point source like a quasar.</p>
<p>High-order interpolation treats the whole sky the same, and has the advantage that it
can be done blindly, with no prior knowledge about which parts of the sky might
be important sources of leakage. This comes at a large performance
penalty, however. And as we have seen, the most accessible such methods, linear and cubic spline
interpolation, do not perfectly eliminate this bias.</p>
<p>If one does know the location
and amplitude of the bright point sources in the sky, however, a much simpler approach
is to simply subtract those point sources in the data before making the map
(<a href=https://arxiv.org/abs/1208.0050>Dunner et al. 2012</a>). This
results in a map free of both the point sources (which can be useful in itself)
and the artifacts they would have sourced. If needed, the source model can then be
projected onto the sky using an uncorrelated noise model $W$ (for example $W=I$)
and added back to the map.
\begin{align}
\hat m_\textrm{srcfree} &= (P^TN^{-1}P)^{-1}P^TN^{-1}(d - \sum_i A_i B_i) \\
\hat m_\textrm{tot} &= \hat m_\textrm{srcfree} + (P^TW^{-1}P)^{-1}P^T W^{-1}\sum_i A_i B_i
\end{align}
Here $A_i$ is the amplitude of source $i$, and $B_i$ is a vector of each sample's
response to an instrument beam centered on the position of that source.
</p>
<h2>3. Inconsistency</h2>
<p>Real telescopes don't have 100% accurate pointing and gain, and even if they did, the
sky is not static and unchanging. These effects mean the telescope effectively sees the
sky jitter around slightly while varying in brightness. There is no room for this in
eq. \ref{eq:model} aside from the noise term $n$, so time-variable instrumental errors
or variable point sources will
lead to exactly the same sort of leakage as
sub-pixel errors do. But since these errors are not sourced by sub-pixel structures,
neither higher-order interpolation nor point source subtraction help here.</p>
<p>The leakage is ultimately caused by two factors:</p>
<ol>
<li>Multiple samples are mapped onto the same pixel, but don't agree on what that
pixel should be.</li>
<li>The correlated noise model causes a non-local response to such errors.</li>
</ol>
<p>This points to two approaches for eliminating the leakage: Model the noise as
uncorrelated, or add more degrees of freedom.</p>
<h3>3.1. White source mapping</h3>
<p>Of course, simply mapping using an uncorrelated noise model would produce a horribly stripy
map in the presence of
realistic correlated noise. However, we don't need everything to be uncorrelated to get
rid of leakage, we just need to decouple the bright point sources from those correlations.
This suggests the following approach:</p>
<ol>
<li>Make a point source free set of samples $d'$ from $d$ by smoothly in-painting the samples in $d$
that hit the bright sources, ideally using maximum likelihood in-painting, but the quality of the
in-painting only matters for the optimality of the method - bad in-painting will not cause a bias.</li>
<li>Make a maximum-likelihood map $\hat m'$ of the in-painted TOD using the full noise model.
\begin{align}
\hat m' &= (P^TN^{-1}P)^{-1}P^T N^{-1} d'
\end{align}
</li>
<li>If the in-painted regions were small compared to the noise correlation length, then
the missing signal $d-d'$ from the previous map will have approximately white noise,
so we can map these bothersome samples using a simple white noise model:
\begin{align}
\Delta m &= (P^TW^{-1}P)^{-1}P^T W^{-1}(d-d')
\end{align}
</li>
<li>Add these maps to get the total, unbiased sky map: $\hat m = \hat m' + \Delta m$.</li>
</ol>
<h3>3.2. Source sub-sampling</h3>
<p>We can eliminate model error by having at least one degree of freedom per sample,
but if we did that to every sample, there would be no averaging down of the noise,
and we wouldn't really get a map, just a reorganized version of the input data. However,
there is nothing wrong with oversampling like this just around strong point sources and
other problematic regions, as long as these regions are small compared to the noise
correlation length. This results in the model
\begin{align}
d &= Pm + Gs + n \\
&= \begin{bmatrix}P&G\end{bmatrix}\begin{bmatrix}m&s\end{bmatrix}^T + n \tag{4} \label{eq:joint}
\end{align}
Here $s$ are the new per-sample degrees of freedom for the samples that hit bright sources,
and $G$ the matrix that maps them 1-to-1 only the corresponding samples. To be explicit,
$a=Gs$ could be implemented as <code>a[mask] = s</code>, and $b=G^Td$ as <code>b = d[mask]</code>,
where mask is True for samples that hit bright sources and False otherwise.</p>
<p>One can then solve jointly for the maximum-likelihood values $\hat m$ and $\hat s$. This can still
be done using equation \ref{eq:mapmaking}, since equation \ref{eq:joint} is still of the same form
as equation \ref{eq:model}. $m$ and $s$ are in isolation not well defined. Some samples are
both covered by pixels from $m$ and by special per-sample degrees of freedom from $m$, so
signal can move freely from one to the other without changing the residuals. This can be avoided by
modifying $P$ to have it skip those samples, but it is not actually necessary, for
the total map
\begin{align}
\hat m_\textrm{tot}&= \hat m + (P^TW^{-1}P)^{-1}P^TW^{-1}\hat s \tag{5} \label{eq:jointsol}
\end{align}
<em>is</em> well-defined in either case, and is the artifact-free map of the sky we are after.</p>
<p>An advantage of this method is that it is very similar to a common way to handle
per-sample cuts in maximum-likelihood map-making. There one also solves for an extra degree
of freedom per cut sample. The only difference is that we here don't throw that part of the
solution away afterwards, but instead project them back into the final map using a white noise model.
cite something here too</p>
<!--
<h4>3.2.1. Why can't this just be applied everywhere?</h4>
<p>At this point one might wonder why one can't just put the whole sky in the mask, and
thus be immune to leakage everywhere. If one were to do this, there would be a degree of
freedom for every sample, so the solution would be free to just match all the raw samples
directly. That is, one would end up with $s = d$, up to degeneracies. The reprojection of
$s$ onto the sky in equation \ref{eq:jointsol} would therefore end up being just a standard
maximum-likelihood map with an uncorrelated noise model, which would be horribly suboptimal
in the presence of realistic noise.</p>
<p>Ultimately it's the constraint that multiple different
samples with multiple different scanning directions map to the same degree of freedom that
allows the correlated noise model to average the noise down. Lifting this constraint in small
regions, like source sub-sampling does, is OK as long as those regions are smaller than the
typical length scale of the correlated noise, but will be progressively less optimal
(more residual correlated noise inside the regions) as they grow beyond that.</p>
-->
<h3>3.3. Source cutting</h3>
<p>Perhaps the conceptually simplest way to avoid model errors from small problematic regions
of the sky is to simply exclude samples hitting those regions from the mapping process. In
practice, however, most correlated noise models are best represented in Fourier space, and fast Fourier
transforms cannot easily handle missing samples. Hence, simply dropping those samples isn't
straightforward. A trick that effectively achieves the same thing is to modify the
pointing matrix to give the bad samples a dummy degree of freedom each instead of pointing
them at the sky (<a href=https://arxiv.org/abs/1208.0050>Dunner et al. 2012</a>),
\begin{align}
d &= \bar P m + Gs + n
\end{align}
Here $G$ has the same definition as in eq. $\ref{eq:joint}$, selecting only the sources hitting
the bright objects, and $\bar P$ is identical to $P$ except that those samples are skipped.
The solution to this is equivalent to $\hat m_\textrm{tot}$ from eq. \ref{eq:jointsol} with
the cut region masked. To see this, notice that replacing $\bar P$ with $P$ will not change
the solution outside the cut area because $\bar P$ already has enough degrees of freedom to
perfectly match the data there, resulting in zero local residual, and adding even more degrees of freedom
in this area will not change that.</p>
<h3>3.4. Summary of the methods</h3>
<p>Table 1 shows a summary of the methods we will investigate. The ideal method would be
blind (no prior knowledge about which parts of the sky have strong signal),
and would handle both sub-pixel errors and inconsistency at low cost, but none of the methods
tick all those boxes, and in particular none of the blind methods can handle inconsistency.</p>
<figure class=wide>
<table class=text>
<tr>
<th></th>
<th>standard</th>
<th>it.lin</th>
<th>it.cubic</th>
<th>bilin</th>
<th>bicubic</th>
<th>srcsub</th>
<th>srcmask</th>
<th>srcwhite</th>
<th>srcsamp</th>
</tr>
<tr>
<th>blind</th>
<td class=good>yes</td>
<td class=good>yes</td>
<td class=good>yes</td>
<td class=good>yes</td>
<td class=good>yes</td>
<td class=bad>no</td>
<td class=bad>no</td>
<td class=bad>no</td>
<td class=bad>no</td>
</tr>
<tr>
<th>subpix</th>
<td class=bad>no</td>
<td class=good>yes</td>
<td class=good>yes</td>
<td class=good>yes</td>
<td class=good>yes</td>
<td class=good>yes</td>
<td class=good>yes</td>
<td class=good>yes</td>
<td class=good>yes</td>
</tr>
<tr>
<th>incon.</th>
<td class=bad>no</td>
<td class=bad>no</td>
<td class=bad>no</td>
<td class=bad>no</td>
<td class=bad>no</td>
<td class=bad>no</td>
<td class=good>yes</td>
<td class=good>yes</td>
<td class=good>yes</td>
</tr>
<tr>
<th>holes</th>
<td class=good>no</th>
<td class=good>no</th>
<td class=good>no</th>
<td class=good>no</th>
<td class=good>no</th>
<td class=good>no</th>
<td class=bad>yes</th>
<td class=good>no</th>
<td class=good>no</th>
</tr>
<tr>
<th>cost</th>
<td class=good>1</td>
<td class=good>1*</td>
<td class=good>1*</td>
<td class=bad>4</td>
<td class=bad>16</td>
<td class=good>1</td>
<td class=good>1</td>
<td class=good>1</td>
<td class=good>1</td>
</tr>
<tr>
<th>compl.</th>
<td class=good>v.low</td>
<td class=medi>high</td>
<td class=bad>v.high</td>
<td class=medi>high</td>
<td class=bad>v.high</td>
<td class=good>low</td>
<td class=good>low</td>
<td class=good>low</td>
<td class=good>low</td>
</tr>
<tr>
<th>pixwin</th>
<td class=neutral>standard</td>
<td class=mixed>small</td>
<td class=mixed>small</td>
<td class=mixed>small</td>
<td class=mixed>small</td>
<td class=neutral>standard</td>
<td class=neutral>standard</td>
<td class=neutral>standard</td>
<td class=neutral>standard</td>
</tr>
<tr>
<th>new</th>
<td class=neutral>no</td>
<td class=mixed>yes</td>
<td class=mixed>yes</td>
<td class=mixed>yes</td>
<td class=mixed>yes</td>
<td class=neutral>no</td>
<td class=neutral>no</td>
<td class=mixed>yes</td>
<td class=mixed>yes</td>
</tr>
</table>
<figcaption>
<dfn>Table 1</dfn>: Summary of the methods discussed in this paper. The rows are: <dfn>blind</dfn>: Whether
the method can be used on the whole map with no prior knowledge of which regions have
large model errors; <dfn>subpix</dfn>: Whether the method targets sub-pixel errors;
<dfn>incon.</dfn>: Whether the method targets inconsistency errors like gain and pointing
errors; <dfn>cost</dfn>: A rough indication of time cost of the method, relative to
standard map-making. Iterative linear/cubic mapmaking naively have a cost factor of
$N_\textrm{iter}$, but in practice this can usually be avoided. <dfn>compl.</dfn>:
A rough indication of how complex the method is to implement (this is somewhat subjective).
<dfn>pixwin</dfn>:
The pixel window the method produces. Higher-order interpolation methods results in
a "smaller" (closer to unity) pixel window. <dfn>new</dfn>: Whether this is a new
method or not. Those marked "yes" are, to my knowledge, not previously published.
</figcaption>
</figure>
<h2>4. Simulations</h2>
<p>I wrote a <a href=https://raw.githubusercontent.com/amaurea/model_error/master/toy.py>minimal maximum-likelihood
mapmaking library</a> in about 300 lines of Python to demonstrate the performance of these
methods, along with <a href=https://raw.githubusercontent.com/amaurea/model_error/master/examples.py>a driver script</a>
that generates the plots in this paper. The library emphasizes simplicity and is neither general
nor optimized for speed, and only depends on <code>numpy</code> and <code>scipy</code> (with
the exception of bilinear and bicubic map-making which also depends on <code>pixell</code>.
That said, the techniques described here have also been implemented in
the full Atacama Cosmology Telescope (ACT) mapmaker, and despite the simplicity of these test cases
the results are representative of running a real map-maker on real data.</p>
<h3>4.1. 2D simulations</h3>
<p>The main simulations use two sets of 243 x 243 equi-spaced samples in a square grid,
the first ordered vertically and the second horizontally (the ordering matters for the
noise correlation structure). These are mapped onto an 81 x 81
pixel grid so that there are 3 x 3 samples per pixel per set. Using this setup I simulated
two signals:</p>
<ol>
<li>A Gaussian point source with a beam standard deviation of 1 pixel and an
amplitude of $2\cdot 10^4$ (500 times higher than the color scale used in the plots in the
results section). This will be representative of e.g. a bright quasar as seen by a high-resolution
CMB telesope like ACT.</li>
<li> A CMB-like field built by simulating a $1/l$ spectrum with an angular
resolution of 1 pixel and an RMS of 1. This will test whether model error
is a problem in the absence of bright sources.</p>
</ol>
<p>I simulated two types of noise model:</p>
<ol>
<li>An uncorrelated noise model $N_\textrm{uncorr} = I$, representing the baseline case
where no leakage is expected.</li>
<li>A correlated noise model $N_\textrm{corr}(f_1,f_2) = \sigma^2 (1 + (f_1/f_\textrm{knee})^\alpha) \delta_{f_1,f_2}$. This is diagonal in Fourier space. $\alpha = -4$ represents an atmosphere-like steeply falling spectrum, meeting a noise floor with RMS $\sigma=1$ at the knee frequency $f_\textrm{knee} = 0.02$ in pixel units.</li>
</ol>
<p>Becuase it is the non-local weighting implied by a correlated noise model that causes leakage,
not the noise itself, one does not actually have to include the noise in the simulations.
To maximize S/N in the figures most of the simulations are noise-free, but I also include some noisy
simulations to show how effectively each method suppresses noise.</p>
<p>I simulated 3 types of model errors:</p>
<ol>
<li>Sub-pixel signal was simulated generating the signals at higher resolution than the
output maps, as described above. This was included in every simulation.</li>
<li>Gain errors were simulated by decreasing/increasing the signal amplitude by
0.5% for the vertical/horizontal data set, representing a total gain mismatch of
1%. This was included in the simulations labeled "gain error".</li>
<li>Pointing errors were simulated by offsetting the horizontal dataset by 0.01
pixels diagonally relative to the vertical dataset. This was includes in the
simulations labeled "pt. error".</li>
</ol>
<h3>4.1. 1D simulations</h3>
<p>I also produced as smaller number of 1-dimensional simulations for Figure 1 and Figure 2.
These used the same noise models and point source signal as the 2-dimensional simulations,
but consisted of just a single scan of length 1100 samples which was mapped onto 100 pixels.</p>
<h2>5. Results</h2>
<p>Figure 3 shows maximum likelihood maps for each method for a strong point source.
In the noise free case (top 3rd of the figure) standard nearest-neighbor mapmaking
with an uncorrelated noise model (first column) results in no leakage as expected, whether in the
presence of normal sub-pixel signal, gain errors or pointing errors. However, in the presence
of realistic correlated noise (mid 3rd of the figure) this method is extremely suboptimal,
resulting in a map completely
dominated by large-scale noise.</p>
<p>Standard map-making with a correlated noise model (second column)
down-weights the noise correctly, but causes an X-shaped pattern of leakage around the
source, regardless of whether the noise is actually present. This only gets worse in
the presence of gain and pointing errors. The size of the leakage relative to
the peak source amplitude is shown in table 2, and is typically of order 3e-4 for
this choice of beam, pixel size and noise model. While
this is small relative to the source itself, it might be much brighter than any
surrounding CMB for a sufficiently bright source.</p>
<p>Full and iterative linear and bicubic interpolation all behave similarly.
In the absence of gain and pointing errors they provide a moderate leakage reduction,
typically by a factor of 3-10, but varying by method and distance from the source.
The computationally heaviest method, bicubic interpolation, provides the best reduction (about a factor
of 20) of these at large distance, but as we saw in figure 2, this comes at the cost of making
things worse near the source.
The iterative approximation to bicubic interpolation might be the best compromise
between low-range and short-range leakage, reducing both to $<0.4 \cdot 10^{-4}$.
However, since these methods only target sub-pixel
errors, it comes as no surprise that they do not help at all with data inconsistencies
like gain and pointing errors, where the X returns in full force.</p>
<p>Source subtraction completely eliminates the artifacts in the ideal case, but
performs no better than the blind interpolation methods in the presence of gain
and pointing errors. The other targeted methods fare better:
Source cutting avoids artifacts even for inconsistent data, at the cost of a hole in the map,
while white source mapping and source sub-sampling avoid both these problems,
eliminating artifacts even in the presence of inconsistent data while still providing
a map of the source.</p>
<figure class=wide>
<table class=pix>
<tr>
<th></th>
<th>standard</th>
<th>standard</th>
<th>it.lin</th>
<th>it.cubic</th>
<th>bilin</th>
<th>bicubic</th>
<th>srcsub</th>
<th>srccut</th>
<th>srcwhite</th>
<th>srcsamp</th>
</tr>
<tr class=sub>
<th></th>
<th>uncorr</th>
<th>corr</th>
<th>corr</th>
<th>corr</th>
<th>corr</th>
<th>corr</th>
<th>corr</th>
<th>corr</th>
<th>corr</th>
<th>corr</th>
</tr>
<tr>
<th>plain</th>
<td><img src=examples/src_2d_uncorr_map_0.png></td>
<td><img src=examples/src_2d_corr_map_0.png></td>
<td><img src=examples/src_2d_corr_itlin_map_0.png></td>
<td><img src=examples/src_2d_corr_itcubic_map_0.png></td>
<td><img src=examples/src_2d_corr_lin_map_0.png></td>
<td><img src=examples/src_2d_corr_cubic_map_0.png></td>
<td><img src=examples/src_2d_corr_srcsub_map_0.png></td>
<td><img src=examples/src_2d_corr_srccut_map_0.png></td>
<td><img src=examples/src_2d_corr_srcwhite_map_0.png></td>
<td><img src=examples/src_2d_corr_srcsamp_map_0.png></td>
</tr>
<tr>
<th>gain error</th>
<td><img src=examples/src_gain_2d_uncorr_map_0.png></td>
<td><img src=examples/src_gain_2d_corr_map_0.png></td>
<td><img src=examples/src_gain_2d_corr_itlin_map_0.png></td>
<td><img src=examples/src_gain_2d_corr_itcubic_map_0.png></td>
<td><img src=examples/src_gain_2d_corr_lin_map_0.png></td>
<td><img src=examples/src_gain_2d_corr_cubic_map_0.png></td>
<td><img src=examples/src_gain_2d_corr_srcsub_map_0.png></td>
<td><img src=examples/src_gain_2d_corr_srccut_map_0.png></td>
<td><img src=examples/src_gain_2d_corr_srcwhite_map_0.png></td>
<td><img src=examples/src_gain_2d_corr_srcsamp_map_0.png></td>
</tr>
<tr>
<th>pt. error</th>
<td><img src=examples/src_ptoff_2d_uncorr_map_0.png></td>
<td><img src=examples/src_ptoff_2d_corr_map_0.png></td>
<td><img src=examples/src_ptoff_2d_corr_itlin_map_0.png></td>
<td><img src=examples/src_ptoff_2d_corr_itcubic_map_0.png></td>
<td><img src=examples/src_ptoff_2d_corr_lin_map_0.png></td>
<td><img src=examples/src_ptoff_2d_corr_cubic_map_0.png></td>
<td><img src=examples/src_ptoff_2d_corr_srcsub_map_0.png></td>
<td><img src=examples/src_ptoff_2d_corr_srccut_map_0.png></td>
<td><img src=examples/src_ptoff_2d_corr_srcwhite_map_0.png></td>
<td><img src=examples/src_ptoff_2d_corr_srcsamp_map_0.png></td>
</tr>
<tr>
<th>plain + noise</th>
<td><img src=examples/src_noise_2d_uncorr_map_0.png></td>
<td><img src=examples/src_noise_2d_corr_map_0.png></td>
<td><img src=examples/src_noise_2d_corr_itlin_map_0.png></td>
<td><img src=examples/src_noise_2d_corr_itcubic_map_0.png></td>
<td><img src=examples/src_noise_2d_corr_lin_map_0.png></td>
<td><img src=examples/src_noise_2d_corr_cubic_map_0.png></td>
<td><img src=examples/src_noise_2d_corr_srcsub_map_0.png></td>
<td><img src=examples/src_noise_2d_corr_srccut_map_0.png></td>
<td><img src=examples/src_noise_2d_corr_srcwhite_map_0.png></td>
<td><img src=examples/src_noise_2d_corr_srcsamp_map_0.png></td>
</tr>
<tr>
<th>gain error + noise</th>
<td><img src=examples/src_noise_gain_2d_uncorr_map_0.png></td>
<td><img src=examples/src_noise_gain_2d_corr_map_0.png></td>
<td><img src=examples/src_noise_gain_2d_corr_itlin_map_0.png></td>
<td><img src=examples/src_noise_gain_2d_corr_itcubic_map_0.png></td>
<td><img src=examples/src_noise_gain_2d_corr_lin_map_0.png></td>
<td><img src=examples/src_noise_gain_2d_corr_cubic_map_0.png></td>
<td><img src=examples/src_noise_gain_2d_corr_srcsub_map_0.png></td>
<td><img src=examples/src_noise_gain_2d_corr_srccut_map_0.png></td>
<td><img src=examples/src_noise_gain_2d_corr_srcwhite_map_0.png></td>
<td><img src=examples/src_noise_gain_2d_corr_srcsamp_map_0.png></td>
</tr>
<tr>
<th>pt. error + noise</th>
<td><img src=examples/src_noise_ptoff_2d_uncorr_map_0.png></td>
<td><img src=examples/src_noise_ptoff_2d_corr_map_0.png></td>
<td><img src=examples/src_noise_ptoff_2d_corr_itlin_map_0.png></td>
<td><img src=examples/src_noise_ptoff_2d_corr_itcubic_map_0.png></td>
<td><img src=examples/src_noise_ptoff_2d_corr_lin_map_0.png></td>
<td><img src=examples/src_noise_ptoff_2d_corr_cubic_map_0.png></td>
<td><img src=examples/src_noise_ptoff_2d_corr_srcsub_map_0.png></td>
<td><img src=examples/src_noise_ptoff_2d_corr_srccut_map_0.png></td>
<td><img src=examples/src_noise_ptoff_2d_corr_srcwhite_map_0.png></td>
<td><img src=examples/src_noise_ptoff_2d_corr_srcsamp_map_0.png></td>
</tr>
<tr>
<th>plain + cmb + noise</th>
<td><img src=examples/src_cmb_noise_2d_uncorr_map_0.png></td>
<td><img src=examples/src_cmb_noise_2d_corr_map_0.png></td>
<td><img src=examples/src_cmb_noise_2d_corr_itlin_map_0.png></td>
<td><img src=examples/src_cmb_noise_2d_corr_itcubic_map_0.png></td>
<td><img src=examples/src_cmb_noise_2d_corr_lin_map_0.png></td>
<td><img src=examples/src_cmb_noise_2d_corr_cubic_map_0.png></td>
<td><img src=examples/src_cmb_noise_2d_corr_srcsub_map_0.png></td>
<td><img src=examples/src_cmb_noise_2d_corr_srccut_map_0.png></td>
<td><img src=examples/src_cmb_noise_2d_corr_srcwhite_map_0.png></td>
<td><img src=examples/src_cmb_noise_2d_corr_srcsamp_map_0.png></td>
</tr>
<tr>
<th>gain error + cmb + noise</th>
<td><img src=examples/src_cmb_noise_gain_2d_uncorr_map_0.png></td>
<td><img src=examples/src_cmb_noise_gain_2d_corr_map_0.png></td>
<td><img src=examples/src_cmb_noise_gain_2d_corr_itlin_map_0.png></td>
<td><img src=examples/src_cmb_noise_gain_2d_corr_itcubic_map_0.png></td>
<td><img src=examples/src_cmb_noise_gain_2d_corr_lin_map_0.png></td>
<td><img src=examples/src_cmb_noise_gain_2d_corr_cubic_map_0.png></td>
<td><img src=examples/src_cmb_noise_gain_2d_corr_srcsub_map_0.png></td>
<td><img src=examples/src_cmb_noise_gain_2d_corr_srccut_map_0.png></td>
<td><img src=examples/src_cmb_noise_gain_2d_corr_srcwhite_map_0.png></td>
<td><img src=examples/src_cmb_noise_gain_2d_corr_srcsamp_map_0.png></td>
</tr>
<tr>
<th>pt. error + cmb + noise</th>
<td><img src=examples/src_cmb_noise_ptoff_2d_uncorr_map_0.png></td>
<td><img src=examples/src_cmb_noise_ptoff_2d_corr_map_0.png></td>
<td><img src=examples/src_cmb_noise_ptoff_2d_corr_itlin_map_0.png></td>
<td><img src=examples/src_cmb_noise_ptoff_2d_corr_itcubic_map_0.png></td>
<td><img src=examples/src_cmb_noise_ptoff_2d_corr_lin_map_0.png></td>
<td><img src=examples/src_cmb_noise_ptoff_2d_corr_cubic_map_0.png></td>
<td><img src=examples/src_cmb_noise_ptoff_2d_corr_srcsub_map_0.png></td>
<td><img src=examples/src_cmb_noise_ptoff_2d_corr_srccut_map_0.png></td>
<td><img src=examples/src_cmb_noise_ptoff_2d_corr_srcwhite_map_0.png></td>
<td><img src=examples/src_cmb_noise_ptoff_2d_corr_srcsamp_map_0.png></td>
</tr>
</table>
<figcaption>
<p><dfn>Figure 3</dfn>: Model errors around a bright point source for the
different methods discussed in this paper, for a simulation with 3x3
samples uniformly distributed inside each pixel for each of a vertical and
horizontal scanning pattern. The top third is for a noise free
simulation, the second third includes strongly correlated noise,
and the bottom third adds a CMB-like smooth signal.
The rows are: <dfn>plain</dfn>: no special conditions; <dfn>gain error</dfn>:
the vertical scanning pattern sees a 1% higher signal than the horizontal
scanning pattern; and <dfn>pt. error</dfn>: the horizontal scanning pattern
sees the signal offset by 0.01 pixel diagonally. The column headers refer to
the different mapmaking algorithms described in the text. The sub-headers
indicate whether a correlated or uncorrelated noise model was used. The
color range is ±$10^{-4}$ of the peak source amplitude.</p>
<!--<p>Standard nearest neighbor mapmaking results in a prominent X-shaped
sub-pixel errors
following the scanning directions as soon as we turn on correlations in the
noise model. With the exception of the bicubic method, the position-independent
methods only modestly reduce these errors. Bicubic mapmaking greatly reduces
the long-distance errors, but introduces some short-distance ringing as the
spline sacrifices fidelity at the edge of the source to match the strong signal
at its center as closely as possible. None of these position-independent methods
help noticeably against non-local bias from gain and pointing errors.</p>-->
</figcaption>
</figure>
<!--
<p>Table 2 gives a more quantitative summary of the amount of leakage from each of these
methods for the noiseless point source case, corresponding to the top 3rd of Figure 3.
The leakage is measured in two regions: near the source (but not inside it), defined as
the area between 5 to 10 pixels from the source center; and far from the source, defined
as the area more than 10 pixels away. The leakage is measured as the highest absolute value
of the map in each region, divided by the peak source amplitude.</p>
<p>For the beam, pixel size and nosie model used here, sub-pixel errors result in a relative
leakage of $\sim 3 \cdot 10^{-4}$. Higher order interpolation in the pointing matrix reduces
this to $0.18 \cdot 10^{-4}$ to $1.79 \cdot 10^{-4}$ depending on the method and distance, with
the iterative bicubic approximation giving relatively good performance
these results.</p>
-->
<figure class=wide>
<table class="text fit">
<tr>
<th></th>
<th>standard</th>
<th>it.lin</th>
<th>it.cubic</th>
<th>bilin</th>
<th>bicubic</th>
<th>srcsub</th>
<th>srccut</th>
<th>srcwhite</th>
<th>srcsamp</th>
</tr>
<tr>
<th>near</th>
<td class=bad >0.866</td>
<td class=medi>0.143</td>
<td class=medi>0.270</td>
<td class=medi>0.172</td>
<td class=bad >1.776</td>
<td class=good>0.000</td>
<td class=good>0.000</td>
<td class=good>0.000</td>
<td class=good>0.000</td>
</tr>
<tr>
<th>far</th>
<td class=bad >3.520</td>
<td class=bad >1.146</td>
<td class=medi>0.392</td>
<td class=bad >0.889</td>
<td class=medi>0.179</td>
<td class=good>0.000</td>
<td class=good>0.000</td>
<td class=good>0.000</td>
<td class=good>0.000</td>
</tr>
<tr>
<th>near + pt.</th>
<td class=bad >2.542</td>
<td class=bad >1.966</td>
<td class=bad >2.057</td>
<td class=bad >2.015</td>
<td class=bad >3.613</td>
<td class=bad >1.910</td>
<td class=good>0.000</td>
<td class=good>0.000</td>
<td class=good>0.000</td>
</tr>
<tr>
<th>far + pt.</th>
<td class=bad >3.605</td>
<td class=bad >1.433</td>
<td class=bad >1.398</td>
<td class=bad >1.439</td>
<td class=bad >1.493</td>
<td class=bad >1.440</td>
<td class=good>0.000</td>
<td class=good>0.000</td>
<td class=good>0.000</td>
</tr>
<tr>
<th>near + gain</th>
<td class=bad >5.665</td>
<td class=bad >4.593</td>
<td class=bad >4.928</td>
<td class=bad >4.605</td>
<td class=bad >6.424</td>
<td class=bad >4.799</td>
<td class=good>0.000</td>
<td class=good>0.000</td>
<td class=good>0.000</td>
</tr>
<tr>
<th>far + gain</th>
<td class=bad >3.990</td>
<td class=bad >3.360</td>
<td class=bad >3.515</td>
<td class=bad >3.350</td>
<td class=bad >3.589</td>
<td class=bad >3.620</td>
<td class=good>0.000</td>
<td class=good>0.000</td>
<td class=good>0.000</td>
</tr>
</table>
<figcaption>
<dfn>Table 2</dfn>: The maximum absolute value of the leakage from a point source
in units of 1e-4 of the peak amplitude, measured in two distance bins: <dfn>near</dfn>:
5-10 pixels from from the center (just outside the point where the beam becomes negligible),
and <dfn>far</dfn>: beyond 10 pixels from the center. For the beam, pixel size and noise model
used here, sub-pixel errors result in O(1e-4) leakage. Higher-order interpolation reduces this
by a factor of 0.5-20 depending on the distance and method (yes, in the worst case it makes it
2x worse near the source), while the targeted methods completely eliminate sub-pixel errors.
Once gain or pointing errors are added, only source cutting, white source mapping and
source sub-sampling provide any meaningful improvement (though this will depend on the
size of the errors chosen).
</figcaption>
</figure>
<p>Figure 4 shows maximum-likelihood maps, residual RMS and deviation from the uncorrelated
solution for each method for a smooth, CMB-like field. It doesn't really make sense to apply
the source-targeted methods when there is no source present, but I still include them with
a dummy point source with zero amplitude in the patch center to see how these methods interact
with the surrounding CMB.</p>
<p>With a smooth, low-dynamic range field like this, the leakage is small
enough that the maps are practically identical. From table 3 and the bottom
3rd of figure 4 we see that the leakage is typically more than 200 times smaller than the
signal itself, even in the presence of gain and pointing errors.</p>
<p>The exception is white
source mapping, which has X-shaped leakage at 3% of the signal strength.
This happens because the
maximum likelihood in-painting of the samples that hit the source is done independently
for the horizontal and vertical data set, which therefore end up disagreeing with
each other in this region. As we have seen before, such disagreement is interpreted
as noise, and spreads out by a noise correlation length. Like the X around strong
point sources, this too is a form of model error bias, but since it's sourced by
the CMB itself rather than a potentially very bright point source, it is never
more than a small fraction of the CMB. A 3%-level contamination in a small fraction
of the sky is too small to matter in most circumstances, but to be safe one can
always use source sub-sampling instead, which doesn't have this issue.
</p>
<figure class=wide>
<table class=pix>
<tr>
<th></th>
<th>standard</th>
<th>standard</th>
<th>it.lin</th>
<th>it.cubic</th>
<th>bilin</th>
<th>bicubic</th>
<th>srcsub</th>
<th>srccut</th>
<th>srcwhite</th>
<th>srcsamp</th>
</tr>
<tr class=sub>
<th></th>
<th>uncorr</th>
<th>corr</th>
<th>corr</th>
<th>corr</th>
<th>corr</th>
<th>corr</th>
<th>corr</th>
<th>corr</th>
<th>corr</th>
<th>corr</th>
</tr>
<tr>
<th>plain</th>
<td><img src=examples/cmb_2d_uncorr_map_0.png></td>
<td><img src=examples/cmb_2d_corr_map_0.png></td>
<td><img src=examples/cmb_2d_corr_itlin_map_0.png></td>
<td><img src=examples/cmb_2d_corr_itcubic_map_0.png></td>
<td><img src=examples/cmb_2d_corr_lin_map_0.png></td>
<td><img src=examples/cmb_2d_corr_cubic_map_0.png></td>
<td><img src=examples/cmb_2d_corr_srcsub_map_0.png></td>
<td><img src=examples/cmb_2d_corr_srccut_map_0.png></td>
<td><img src=examples/cmb_2d_corr_srcwhite_map_0.png></td>
<td><img src=examples/cmb_2d_corr_srcsamp_map_0.png></td>
</tr>
<tr>
<th>gain error</th>
<td><img src=examples/cmb_gain_2d_uncorr_map_0.png></td>
<td><img src=examples/cmb_gain_2d_corr_map_0.png></td>
<td><img src=examples/cmb_gain_2d_corr_itlin_map_0.png></td>
<td><img src=examples/cmb_gain_2d_corr_itcubic_map_0.png></td>
<td><img src=examples/cmb_gain_2d_corr_lin_map_0.png></td>
<td><img src=examples/cmb_gain_2d_corr_cubic_map_0.png></td>
<td><img src=examples/cmb_gain_2d_corr_srcsub_map_0.png></td>
<td><img src=examples/cmb_gain_2d_corr_srccut_map_0.png></td>
<td><img src=examples/cmb_gain_2d_corr_srcwhite_map_0.png></td>
<td><img src=examples/cmb_gain_2d_corr_srcsamp_map_0.png></td>
</tr>
<tr>
<th>pt. error</th>
<td><img src=examples/cmb_ptoff_2d_uncorr_map_0.png></td>
<td><img src=examples/cmb_ptoff_2d_corr_map_0.png></td>
<td><img src=examples/cmb_ptoff_2d_corr_itlin_map_0.png></td>
<td><img src=examples/cmb_ptoff_2d_corr_itcubic_map_0.png></td>
<td><img src=examples/cmb_ptoff_2d_corr_lin_map_0.png></td>
<td><img src=examples/cmb_ptoff_2d_corr_cubic_map_0.png></td>
<td><img src=examples/cmb_ptoff_2d_corr_srcsub_map_0.png></td>
<td><img src=examples/cmb_ptoff_2d_corr_srccut_map_0.png></td>
<td><img src=examples/cmb_ptoff_2d_corr_srcwhite_map_0.png></td>
<td><img src=examples/cmb_ptoff_2d_corr_srcsamp_map_0.png></td>
</tr>
<tr>
<th>plain resid</th>
<td><img src=examples/cmb_2d_uncorr_map_1.png></td>
<td><img src=examples/cmb_2d_corr_map_1.png></td>
<td><img src=examples/cmb_2d_corr_itlin_map_1.png></td>
<td><img src=examples/cmb_2d_corr_itcubic_map_1.png></td>
<td><img src=examples/cmb_2d_corr_lin_map_1.png></td>
<td><img src=examples/cmb_2d_corr_cubic_map_1.png></td>
<td><img src=examples/cmb_2d_corr_srcsub_map_1.png></td>
<td><img src=examples/cmb_2d_corr_srccut_map_1.png></td>
<td><img src=examples/cmb_2d_corr_srcwhite_map_1.png></td>