Biased Cls when using weighted maps

[tilmantroester]

I have an issue where namaster seems to fail to return unbiased estimates of the power spectrum if a weight map is being used as the mask. I don't know if that's a fundamental issue (doubtful, since support for weight maps is stated in a couple of places) or I messed up somewhere (entirely possible but the same code works for binary masks). Apodisation doesn't seem to help.

The setup is as follows (full code at the end):

• The weight map has pixel weights drawn uniformly between 0 and 1, and then restricted to a circular footprint of fsky = 0.1
• The simulated maps are full-sky noise maps for T, Q, U, with pixel variance \sigma^2_X, X=T,Q,U, \sigma^2_T = 1.5, \sigma^2_Q=\sigma^2_U = 1
• Coupled and decoupled Cls are computed using the workspace setup using the weight map.
• The expected Cls are computed using nmt_workspace.decouple_cell(nmt_workspace.couple_cell(Cl_theory))

For full-sky (top panel) and binary mask (bottom panel, just the circular footprint), this works as expected:

When using the weight map, the measured, decoupled Cls don't agree with the analytic values anymore:

Code:

import numpy as np
import healpy
import pymaster as nmt

import matplotlib.pyplot as plt

nside = 256
n_pix = healpy.nside2npix(nside)

pixel_var_T = 1.5
pixel_var_Q = pixel_var_U = 1.0

w = np.random.rand(n_pix)

fsky = 0.1
mask = np.ones(n_pix)
mask[int(fsky*n_pix):] = 0

mask *= w

# mask = nmt.mask_apodization(mask,
#                             0.5, apotype="Smooth")

f_0 = nmt.NmtField(mask, None, spin=0)
f_2 = nmt.NmtField(mask, None, spin=2)

delta_ell = 10
nmt_bins = nmt.NmtBin.from_nside_linear(nside=nside,
                                        nlb=delta_ell)
ell_eff = nmt_bins.get_effective_ells()

nmt_00_workspace = nmt.NmtWorkspace()
nmt_00_workspace.compute_coupling_matrix(fl1=f_0,
                                         fl2=f_0,
                                         bins=nmt_bins)

nmt_22_workspace = nmt.NmtWorkspace()
nmt_22_workspace.compute_coupling_matrix(fl1=f_2,
                                         fl2=f_2,
                                         bins=nmt_bins)

def simulate_map():
    T = np.random.normal(size=n_pix)*np.sqrt(pixel_var_T)
    Q = np.random.normal(size=n_pix)*np.sqrt(pixel_var_Q)
    U = np.random.normal(size=n_pix)*np.sqrt(pixel_var_U)

    f_rand_0 = nmt.NmtField(mask, [T])
    f_rand_2 = nmt.NmtField(mask, [Q, U])

    Cls_fullsky = healpy.anafast([T,Q,U], nspec=2)

    Cls_00_coupled = nmt.compute_coupled_cell(f_rand_0, f_rand_0)
    Cls_22_coupled = nmt.compute_coupled_cell(f_rand_2, f_rand_2)

    Cls_00_decoupled = nmt_00_workspace.decouple_cell(Cls_00_coupled)
    Cls_22_decoupled = nmt_22_workspace.decouple_cell(Cls_22_coupled)

    return Cls_fullsky[[0,1]], Cls_00_decoupled[0], Cls_22_decoupled[0]

n_ell = 3*nside
n_sim = 10

TT = np.empty((n_sim, n_ell))
EE = np.empty((n_sim, n_ell))

TT_decoupled = np.empty((n_sim, len(ell_eff)))
EE_decoupled = np.empty((n_sim, len(ell_eff)))

for i in tqdm.trange(n_sim):
    (TT[i], EE[i]), TT_decoupled[i], EE_decoupled[i] = simulate_map()

TT_mean = np.mean(TT, axis=0)
TT_err = np.sqrt(np.var(TT, axis=0, ddof=1)/n_sim)

EE_mean = np.mean(EE, axis=0)
EE_err = np.sqrt(np.var(EE, axis=0, ddof=1)/n_sim)

TT_decoupled_mean = np.mean(TT_decoupled, axis=0)
TT_decoupled_err = np.sqrt(np.var(TT_decoupled, axis=0, ddof=1)/n_sim)

EE_decoupled_mean = np.mean(EE_decoupled, axis=0)
EE_decoupled_err = np.sqrt(np.var(EE_decoupled, axis=0, ddof=1)/n_sim)

Cl_TT = pixel_var_T * healpy.nside2pixarea(nside) * np.ones_like(ell)
Cl_EE = (pixel_var_Q + pixel_var_U)/2 * healpy.nside2pixarea(nside) * np.ones_like(ell)

Cl_TT_decoupled = nmt_00_workspace.decouple_cell(nmt_00_workspace.couple_cell([Cl_TT]))
Cl_EE_decoupled = nmt_22_workspace.decouple_cell(nmt_22_workspace.couple_cell([Cl_EE, 0*Cl_EE, 0*Cl_EE, Cl_EE]))

fig, ax = plt.subplots(2, 1, sharex=True, figsize=(5, 4))
fig.subplots_adjust(hspace=0)

ax[0].errorbar(ell_eff, 
               nmt_bins.bin_cell(TT_mean),
               nmt_bins.bin_cell(TT_err)/np.sqrt(delta_ell),
               c="C0", ls="none", marker=".", label="Fullsky anafast TT")
ax[0].errorbar(ell_eff,
               nmt_bins.bin_cell(EE_mean),
               nmt_bins.bin_cell(EE_err)/np.sqrt(delta_ell),
               c="C1", ls="none", marker=".", label="Fullsky anafast EE")
ax[0].plot(ell_eff, nmt_bins.bin_cell(Cl_TT),
           c="C0", label=r"$\sigma_T^2 A_{\rm pix}$")
ax[0].plot(ell_eff, nmt_bins.bin_cell(Cl_EE),
           c="C1", label=r"$(\sigma_Q^2+\sigma_U^2)/2 A_{\rm pix}$")

ax[1].plot(ell_eff, Cl_TT_decoupled[0],
           c="C0", ls="-", label="Namaster, expected TT")
ax[1].errorbar(ell_eff, TT_decoupled_mean, TT_decoupled_err,
               c="C0", ls="none", marker=".", label="Namaster, measured TT")

ax[1].plot(ell_eff, Cl_EE_decoupled[0],
           c="C1", ls="-", label="Namaster, expected EE")
ax[1].errorbar(ell_eff, EE_decoupled_mean, EE_decoupled_err,
               c="C1", ls="none", marker=".",  label="Namaster, measured EE")

ax[0].legend(ncol=2, fontsize="small")
ax[1].legend(ncol=2, fontsize="small")

[a.set_ylim(top=4e-5) for a in ax]
[a.set_ylabel(r"$C_\ell$") for a in ax]
ax[1].set_xlabel(r"$\ell$")

fig.dpi = 150
fig.suptitle("Weights")
fig.savefig("plots/namaster_noise_test_weights.png")

[damonge]

@tilmantroester in this case the mask (and the fields) you're using are complete non-band-limited (i.e. the power spectrum of the mask is constant at high ell), so the fact that namaster only computes mode-coupling up to a given ell_max coupled with the numerical errors of healpix means that the resulting MCM is probably inaccurate.

Maybe try using a nicer mask (e.g. a smoothed version of this one)?

[xzackli]

@damonge is correct here that the mask is pathological, and reducing the power at small scales will make mode-coupling accurate again. This is the usual signal processing situation: input signals which are not bandlimited suffer from aliasing, so you have to apply a low-pass filter.

Here's a row of the mode-coupling matrix comparing a non-band-limited mask, and a nice smoothed mask, plotted in log-log. This describes the contributions to the pseudo-spectrum at ell=200 from the underlying spectrum. The mode-coupling matrix elements for the white noise mask are increasing with ell, because the power spectrum of the mask is constant.

Simply scaling the harmonic coefficients of the mask by sqrt(ell) past some high multipole is sufficient for me to recover unbiased spectra again.

[tilmantroester]

The issue persists when using a real weight map. It's a shear weight map, so it's not clear to me how that can be smoothed while conserving the shear weight information appropriately. Not using nmt_workspace.couple_cell but instead scaling the theoretical noise spectrum by the mean of the square of the weight map solved the problem, however:

This can be closed I think. A heuristic that checks if a mask is problematic and warns the user might be useful though.

[damonge]

Yes, this is due to the finite ell_max up to which the MCM is computed, whereas the input power spectrum has power up to ell=infinity in principle. Multiplying by <mask^2> here is the exact effect of a mask on a flat power spectrum, so that's bound to do better (this is incidentally what we did in https://arxiv.org/abs/2010.09717).

In practice we have never seen this be a problem for shear-like signal power spectra, where the signal decreases with scale relatively fast, but this is something we may need to test again for LSST when the error bars are smaller.

I'll close this for now then, but please feel free to reopen in any related issues appear.