s2let_hpxtest.m with too large error

[flying-gwx]

Hello, I want to use s2let in Healpix sampling scheme. I build the s2let with Healpix_3.82(for lib).

Here is the code in s2let_hpxtest.m I used:

clear all;
close all;

setenv('HEALPIX','/home/gaowenxuan/Code/Healpix_3.82')
setenv('HEALPIXDATA','/home/gaowenxuan/Code/Healpix_3.82/data')
% Main parameters
L = 64;
nside = 128;
B = 2;
J_min = 4;
J = s2let_jmax(L, B);
npix = 12*nside*nside;
Spin = 0;

disp('Generates band-limited function')
flm = zeros(L^2,1); 
flm = rand(size(flm)) + sqrt(-1)*rand(size(flm));
flm = 2.*(flm - (1+sqrt(-1))./2);

disp('Constraint on flms to generate real signal')
for el = 0:L-1
   ind = el*el + el + 1;
   flm(ind,1) = real(flm(ind,1));
   for m = 1:el
      ind_pm = el*el + el + m + 1;
      ind_nm = el*el + el - m + 1;
      flm(ind_nm,1) = (-1)^m * conj(flm(ind_pm,1));
   end  
end


disp('Perform spherical harmonic decomposition with default parameters')
f = s2let_hpx_alm2map(flm, nside, 'L', L);
flm_rec = s2let_hpx_map2alm(f, 'L', L);
default = max(max(abs(flm-flm_rec)))


disp('Perform spherical harmonic decomposition with custom parameters')
f = s2let_hpx_alm2map(flm, nside, 'L', L);
flm_rec = s2let_hpx_map2alm(f, 'nside', nside, 'L', L);
custom = max(max(abs(flm-flm_rec)))
 

disp('Perform axisym wavelet transform with default parameters')
[f_wav, f_scal] = s2let_transform_axisym_analysis_hpx(f);
f_rec = s2let_transform_axisym_synthesis_hpx(f_wav, f_scal);
default = max(max(abs(f-f_rec)))


disp('Perform axisym wavelet transform with custom parameters')
[f_wav, f_scal] = s2let_transform_axisym_analysis_hpx(f, 'B', B, 'L', L, 'J_min', J_min);
f_rec = s2let_transform_axisym_synthesis_hpx(f_wav, f_scal, 'B', B, 'L', L, 'J_min', J_min);
custom = max(max(abs(f-f_rec)))

Here is the output:

Generates band-limited function
Constraint on flms to generate real signal
Perform spherical harmonic decomposition with default parameters

default =

    4.6795

Perform spherical harmonic decomposition with custom parameters

custom =

    4.6795

Perform axisym wavelet transform with default parameters

default =

  963.4893

Perform axisym wavelet transform with custom parameters

custom =

  953.6840

I read the s2let paper, the error sholud be less than 10^-6 for max(|flm - flm_rec|). My error is too large(4.6795).

1. In paper S2LET: A code to perform fast wavelet analysis on the sphere, there exists iteration, but I didn't find it in matlab code s2let_hpx_alm2map or s2let_transform_axisym_analysis_hpx, is there any thing I missed?
2. I found that the Heaplix version is different with the original code (healpix2.20a), does this really matter?

[jasonmcewen]

The forward and inverse wavelet transform, when starting with spherical harmonic coefficients and going back to spherical harmonic coefficients, should indeed reach numerical precision. However, the spherical harmonic transforms with Healpix are highly approximate, so you'll see a large overall error if comparing signals before and after a wavelet transform in pixel space. One way to improve Healpix accuracy is to iterate the Healpix spherical harmonic transforms.

Note also that we're in the process of re-writing ssht, so3 and s2let in JAX to support GPU acceleration and automatic differentiation. This is still a work in progress but the corresponding s2fft and s2wav codes are already public.

[flying-gwx]

Thanks for your reply. My statement may not be clear and let me explain my problem.

The forward and inverse spherical harmonic transform with Healpix, which starting with spherical harmonic coefficients and going back to spherical harmonic coefficients, do not reach numerical precision in spherical harmonic space. The error of 4.6795 is actually in spherical harmonic space and it is too large comparing with the paper(alound 10-6 in spherical harmonic space).

I totally understand the Healpix can result in error comparing signals before and after a wavelet transform in pixel space.
However, the error is too large even in spherical harmonic space.

[jasonmcewen]

Apologies @flying-gwx , we're not able to provide further support for the legacy codes. I would recommend using s2fft and s2wav. Both can support HEALPix, as well as equiangular sampling schemes.