Weights vs Interations

[LeeoBianchi]

Hi @xzackli,

@ziotom78 pointed you as the right person to ask this :)

While working for my thesis on the MPI support for Healpix.jl I faced the problem of how to deal with the iterations in the map2alm transform. My idea was to implement the usage of weights to get rid of iterative routine.

Moreover, while digging into the libsharp2 functions and classes I've noticed that it should be possible to add ring weights to the sharp_geom_info while building it, instead of applying the weights directly to the map as it is implemented now, which sounded like a much more natural solution to me. Am I missing any counter-indication to do this?

Moreover, which is the most correct way to test an sht's accuracy? At the moment I've just tried to transform a weighted map to the harmonic space and then back, computing the variance wrt the original map. But I can only notice a VERY slight difference between a transform with let's say 5 iteration and one with no iterations and no weights and one with weights and no iterations.

Thanks for your help, Leo

[xzackli]

Hey @LeeoBianchi,

There is already some basic support for pixel weights, which you may have seen. These are better than ring weights and involve solving some very large system for weights on each pixel rather than each ring. These "full pixel weights" recover functions with band limit $\ell_{max} \sim (3/2) \times \texttt{nside}$ to machine precision; see this github isuse by @mreineck. The interface for these weights could be improved in Healpix.jl, since it could automatically download weights as the pixwin function does.

A typical way to assess accuracy is the maximum absolute error between the true and transformed (map2alm and back) map. It'd be interesting to see a plot, though I imagine Martin has one somewhere. I believe the map2alm unit tests in Healpix.jl will fail without the proper number of iterations if you want a test on the resulting $a_{\ell m}$ coefficients.

[mreineck]

Pixel weights are great, but there is an important caveat: they ONLY can work properly if your original map is strictly band limited at lmax=1.5*Nside. Realistic maps typically never are band limited, so I'd rather recommend the more dynamic approach implemented in https://github.com/healpy/healpy/pull/734. It stops as soon as the desired accuracy is reached, and if it is not reached, it gives feedback on how good/bad the answer is.

To be honest, I cannot recommend map analysis with weights any more since this new approach exists... your results may be really bad (depending on how well-behaved your input map is), and you are not warned about it.

[xzackli]

The dynamic approach sounds like a great thing to implement. It looks like Martin has provided some pseudocode within the link, and there is a reference implementation in healpy to test against.

[mreineck]

It's also already implemented in ducc0, so I can add it to the wrapped functions :-)

[LeeoBianchi]

Hi again @xzackli @mreineck, Thank you very much for your replies.

The dynamic approach sounds great indeed. I am only scared of it in the MPI parallel setting, since to me it looks like is requiring some extra MPI communication, at least to compute the a_lm norm to compare with the tolerance.

Maybe a user-fixed iteration number could be better in this specific case? Since it avoids the MPI communication for the dot product while each subset of the (residual) map can be subtracted locally on each MPI task (?)

[mreineck]

Yes, the dynamic approach and MPI are quite difficult to bring together ... not something I'd recommend as a side project :-/

I guess the best solution depends strongly on the following parameters:

• how large are your targetted Nside and lmax?
• are your input maps band limited?
• how accurate (in terms of L2 error) do you need the result to be?
• how many cores do your compute nodes have?
• do you have the option of analyzing several maps simultaneously?

I'd say that in the overwhelming majority of cases simple multi-threaded SHTs without MPI will be the best option (and as far as I have heard, this is also planned for future Commander versions).

One more small comment:

Moreover, which is the most correct way to test an sht's accuracy? At the moment I've just tried to transform a weighted map to the harmonic space and then back, computing the variance wrt the original map. But I can only notice a VERY slight difference between a transform with let's say 5 iteration and one with no iterations and no weights and one with weights and no iterations.

Then you have probably chosen a very "nice" input map with steeply falling power spectrum. If you add some amount of white noise to the map and repeat the experiment, you will notice very clear differences between input and output maps, which will probably depend on the number of iterations, but will never go away, no matter how many iterations you use.

[LeeoBianchi]

Yes, the dynamic approach and MPI are quite difficult to bring together ... not something I'd recommend as a side project :-/

I agree, I'll just keep it as a last "if-time" optional for my MPI implementation, and maybe implement pixel weights in the current Healpix.jl fashion instead if really needed.

I guess the best solution depends strongly on the following parameters:

* how large are your targetted `Nside` and `lmax`?

* are your input maps band limited?

* how accurate (in terms of L2 error) do you need the result to be?

* how many cores do your compute nodes have?

* do you have the option of analyzing several maps simultaneously?

This is very useful and I'll have to come back to this when facing a more real case later on in my project. Thanks!

Then you have probably chosen a very "nice" input map with steeply falling power spectrum. If you add some amount of white noise to the map and repeat the experiment, you will notice very clear differences between input and output maps, which will probably depend on the number of iterations, but will never go away, no matter how many iterations you use.

I see! This is it. In fact with fully white noise maps I get a quite high rms between original and output map that does not go away with weights nor high iteration number!

[mreineck]

I see! This is it. In fact with fully white noise maps I get a quite high rms between original and output map that does not go away with weights nor high iteration number!

That's an unfortunate property of spherical geometry. In Cartesian space you can FFT a regular grid containing white noise, and after an inverse FFT you get the input array back (up to very small numerical errors). On the sphere, no grid exists where this can be done ... only the other direction is possible. In other words, there exist grids for which map2alm(alm2map(input)) will always match input, no matter what the values in input were, but there are no grids where you can guarantee alm2map(map2alm(input)) == input.

[LeeoBianchi]

I see! This is it. In fact with fully white noise maps I get a quite high rms between original and output map that does not go away with weights nor high iteration number!

That's an unfortunate property of spherical geometry. In Cartesian space you can FFT a regular grid containing white noise, and after an inverse FFT you get the input array back (up to very small numerical errors). On the sphere, no grid exists where this can be done ... only the other direction is possible. In other words, there exist grids for which map2alm(alm2map(input)) will always match input, no matter what the values in input were, but there are no grids where you can guarantee alm2map(map2alm(input)) == input.

Sure. And If I'm not wrong the same applies to the corresponding adjoint operators, but with a swapped role for map and alm arguments i.e.: I cannot guarantee that adjoint_alm2map(adjoint_map2alm(input_alm)) == input_alm while there exist grids where the other way around is guaranteed to hold true, right?

[mreineck]

I'm not really sure if this is correct, I never thought about this much. But if we assume a combination of lmax and grid geometry where map2alm(alm2map(input_alm)) will always match input_alm, then adjoint_map2alm(adjoint_alm2map(input_map)) cannot generally match input_map, since the map vectors will be longer than the alm vectors,and the adjoint_alm2map operation will destroy information.