allow lmax greater than the exact bandlimit

[zatkins2]

Hi @mreineck,

Just wanted to follow-up in a dedicated place to our discussion here. I guess I have two questions, one more specific and one more general.

More specific

Perhaps naively, I thought the maximum lmax of SHT analysis was 1 more than what's calculated here. For instance, for a regular rectangular pixelization (e.g. CC or F1) with 0.5' pixel size (used in ACT), the ducc lmax is 21599, whereas the Nyquist frequency of such a pixelization is 21600 (with k_nyq = 0.5 * 2*pi / resolution). Perhaps this is related to the novel algorithm you implement (or perhaps my intuition is just wrong)?

More general

Again, unless there are new considerations from your algorithm that go beyond a ``classic" quadrature rule, I'm not sure I agree with your statement that

all you will see at these higher moments is aliasing from moments <=lmax and no real information about what's "really" happening at the high frequencies, because the chosen grid simply cannot resolve them

Considering a regular rectangular pixelization, don't rings farther from the equator sample frequencies beyond lmax in the x-direction? Not sure they could be usefully extracted with an SHT as I agree there will be more and more aliasing as you go farther beyond lmax (really my silly interest is to go 1 beyond lmax, as above, if nothing else because round numbers are nice 😛).

[mreineck]

Hi @zatkins2!

Concerning the more specific question: This is a bit tricky to answer because I guess there are two different concepts of lmax in play. Let's assume a Gauss-Legendre grid for a moment: if you want to store information up to (and including) lmax in it, you need lmax+1 rings. You can still analyze the grid with lmax+1, and you will faithfully get the moment at lmax+1 back - which is (by construction) zero. I'm pretty sure that the same will not work for arbitrary moments at lmax+1 if you have a GL grid with lmax+1 rings.

Concerning the more general question, I have put together a very simple demo script which generates a GL map from alm up to a desired lmax and then tries to analyze it for moments above this limit. Following your argument, these higher moments should be close to zero (we didn't put anything in at these moments after all!), but instead they are perfect mirrors of the lower moments.

I admit that I haven't tried other grids yet; I'll have to make a few adjustments to get this to work, but I'll post my results here!

import ducc0
import numpy as np
from time import time

rng = np.random.default_rng(48)

# helper functions

def nalm(lmax, mmax):
    return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)

def random_alm(lmax, mmax, spin, ncomp):
    res = rng.uniform(-1., 1., (ncomp, nalm(lmax, mmax))) \
     + 1j*rng.uniform(-1., 1., (ncomp, nalm(lmax, mmax)))
    # make a_lm with m==0 real-valued
    res[:, 0:lmax+1].imag = 0.
    ofs=0
    for s in range(spin):
        res[:, ofs:ofs+spin-s] = 0.
        ofs += lmax+1-s
    return res

# test function:
# - create a_lm with lmax=lmax, mmax=0, 
# - put them on a suitable GL grid
# - analyze the grid using analysis_2d, resulting in alm2
# - brute-force "analyze" the grid (up to lmax2) using GL weights and
#   adjoint_synthesis_2d, resulting in alm3
def test(lmax, lmax2):
    geometry= "GL"
    nrings = lmax+1

    # randoom a_lm
    alm = random_alm(lmax, 0, 0, 1)

    # go to minimum GL map
    map = ducc0.sht.experimental.synthesis_2d(alm=alm, lmax=lmax, mmax=0, spin=0, geometry=geometry,ntheta=nrings, nphi=1)

    # standard analysis up to (and including) lmax
    alm2 = ducc0.sht.experimental.analysis_2d(map=map, lmax=lmax, mmax=0, spin=0, geometry=geometry)

    # apply GL weights
    wgt = ducc0.sht.experimental.get_gridweights(geometry, nrings)
    mapx = map*wgt.reshape((1,-1,1))

    # non-standard analysis up to (and including) lmax2
    alm3 = ducc0.sht.experimental.adjoint_synthesis_2d(lmax=lmax2, mmax=0, spin=0, map=mapx, nthreads=1, geometry=geometry)

    # plot the results
    import matplotlib.pyplot as plt
    plt.plot(np.abs(alm2[0]),label="analysis_2d")
    plt.plot(np.abs(alm3[0]),linestyle="dashed", label="nonstandard analysis")
    plt.plot(np.abs(alm[0]),linestyle="dotted", label="input")
    plt.legend()
    plt.show()

test(lmax=10, lmax2=20)

[mreineck]

Here is a slightly different example which also works on non-GL grids. The idea here is to start with a grid that can (presumably) hold moments <= lmax_grid, then synthesize a_lm at a higher moment onto that grid,and finally look at the analyzed a_lm. Optimally, they should be identically zero.

import ducc0
import numpy as np

rng = np.random.default_rng(48)

# helper functions

def nalm(lmax, mmax):
    return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)

def random_alm(lmax, mmax, spin, ncomp):
    res = rng.uniform(-1., 1., (ncomp, nalm(lmax, mmax))) \
     + 1j*rng.uniform(-1., 1., (ncomp, nalm(lmax, mmax)))
    # make a_lm with m==0 real-valued
    res[:, 0:lmax+1].imag = 0.
    ofs=0
    for s in range(spin):
        res[:, ofs:ofs+spin-s] = 0.
        ofs += lmax+1-s
    return res

# test function:
# - create a_lm wwhich are zero except at lmax=lmax, mmax=0, 
# - put them on a grid with a number of rings required for a band limit of lmax2
# - analyze the grid using analysis_2d
# - show resulting a_lm
def test2(geometry, lmax_grid, lmax_data):

    nrings = lmax_grid+1
    if geometry=="CC":
        nrings = lmax_grid+2
    elif geometry=="DH":
        nrings = 2*lmax_grid+2
    elif geometry=="F2":
        nrings = 2*lmax_grid+1

    # random a_lm at (lmax_data,0)
    alm = random_alm(lmax_data, 0, 0, 1)
    alm *= 0
    alm[0,lmax_data] = 1

    # go to minimum map
    map = ducc0.sht.experimental.synthesis_2d(alm=alm, lmax=lmax_data, mmax=0, spin=0, geometry=geometry,ntheta=nrings, nphi=1)

    # standard analysis up to (and including) lmax_grid
    alm2 = ducc0.sht.experimental.analysis_2d(map=map, lmax=lmax_grid, mmax=0, spin=0, geometry=geometry)

    # plot the results
    import matplotlib.pyplot as plt
    plt.plot(np.abs(alm2[0]), label="result")
    plt.legend()
    plt.show()

test2("CC", lmax_data=12, lmax_grid=10)

As you can see from the example, the "too high" moments cause severe aliasing; in other words, the spherical analysis interprets moments beyond the Nyquist frequency as other moments below Nyquist. This effect is bad, and it won't go away if we allow analyzing maps to a "higher than officially allowed" lmax.

Please note that using lmax_data==lmax_grid+1 seems to work OK (at least in the cases I tried). But I haven't convinced myself yet that this is generally true...