Issue in the Conjugate Gradient algorithm (destriper) #322
Description
Using litebird_sim internal destriper, I noticed that the CG algorithm is behaving weirdly when playing with values of the threshold (i.e. the minimum discrepancy between the estimated baselines and the baselines deduced by the destriped map) lower than the default of 1e-7.
threshold = 1e-8 # default is 1e-7
destriper_params = lbs.DestriperParameters(
output_coordinate_system=lbs.coordinates.CoordinateSystem.Galactic,
samples_per_baseline=samples_per_baseline,
iter_max=100,
threshold=threshold,
use_preconditioner=True,
)
Simulation features:
- channel M1-140
- 2 detectors
- 1 year of observation
- white noise level as in the post ptep sims net_ukrts=18.29
- 1/f_knee freq=100mHz
I tested CMB + noise, serially and using various MPI tasks, with different noise realizations (different seeds), and threshold values. I observed:
- for the first 5 iterations, CG effectively minimized the residuals to ~7e-8
- from the 6th iteration onward, CG started oscillating, with residuals increasing (e.g., for
iter_max=100, I gotDestriper CG iteration 100/100, stopping factor: 2.306e+02).
However the issue is "contained" in the sense that the baselines are updated only for lower values of the residuals, therefore baselines are updated up to the 5th iteration but never later. Here in destriper.py:
cur_stopping_factor = _get_stopping_factor(new_r)
if (not best_stopping_factor) or cur_stopping_factor < best_stopping_factor:
best_stopping_factor = cur_stopping_factor
for cur_best_x_k, x_k in zip(best_x, x):
cur_best_x_k[:] = x_k
This means that setting the threshold value below the default of 1e-7 is producing the exact same result as if using threshold=1e-7, but possibly taking more time (depending on the number of iterations set, iter_max=100).