[feature] Lomb-Scargle Periodograms for Multiband Data

[JulienPeloton]

Hello,

I’m amazed by the performance of the fastnifty method for Lomb-Scargle analysis, and I have started to use it everywhere for LSST analysis (with millions of periods to extract)! I was wondering if there are any plans to add this feature for the multiband version?

[lgarrison]

I'm glad you've found the package useful! For multi-band periodograms, I think it would be great to support, although I haven't looked into it before now. Skimming the Astropy docs & source code, it seems like there's a "fast" implementation and a "flexible" one. The fast one seems like it might "just" be a wrapper around single-band Lomb-Scargle (https://github.com/astropy/astropy/blob/main/astropy/timeseries/periodograms/lombscargle_multiband/implementations/mbfast_impl.py). In other words, you might be able to use the Astropy interface with sb_method="fastnifty" and get a speedup.

The code for the flexible method is more complicated! https://github.com/astropy/astropy/blob/main/astropy/timeseries/periodograms/lombscargle_multiband/implementations/mbflex_impl.py . I can't tell immediately if a NUFFT would help there.

Do you know if you need the fast or flexible method? If fast is acceptable, it would be good to know if the fastnifty backend works for you.

Finally, I'll note that an enhancement that would probably help for millions of periodograms is the ability to thread over multiple periodograms. Right now, we only support threading over periodograms with the same observation times, but we should be able to support different observation times, too. Multiprocessing can do this, but probably not as well as threads.

[JulienPeloton]

Thanks for the detailed answer!

In other words, you might be able to use the Astropy interface with sb_method="fastnifty" and get a speedup.

Indeed, good catch. It works out of the box with the default configuration, with a ~5x speed-up compare to astropy implementation.

Do you know if you need the fast or flexible method? If fast is acceptable, it would be good to know if the fastnifty backend works for you.

Fast is acceptable so far, the results I get are completely decent. However playing a bit with the parameters, I found that setting sb_method="fastnifty" does not allow to have the number of frequency terms (common and per-band) higher than 1:

model = LombScargleMultiband(
    time, mag, filters, error, nterms_base=2, nterms_band=1
)

frequency, power = model.autopower(method="fast", sb_method="fastnifty")
---------------------------------------------------------------------------   
ValueError: nterms != 1 only supported with 'chi2' or 'fastchi2' methods

https://github.com/astropy/astropy/blob/bf1103e5061d4a13c928afcaa66865154f5d9afe/astropy/timeseries/periodograms/lombscargle/implementations/main.py#L221-L226

but definitely not a showstopper at this stage of my work (I rather want speed than precision).

Finally, I'll note that an enhancement that would probably help for millions of periodograms is the ability to thread over multiple periodograms. Right now, we only support threading over periodograms with the same observation times, but we should be able to support different observation times, too. Multiprocessing can do this, but probably not as well as threads.

Currently I use naive distributed processing to parallelize the load (Spark) -- so each periodogram is treated independently from the others. But I agree threading over periodograms with different observation times would be great (as probably many steps are re-used internally) :-)

[lgarrison]

Great! I'm glad that's a partial solution, anyway. I'll keep this issue open to remind us about the flexible method. And I'm glad the package is playing nice with Spark.

Even single-band nifty-ls doesn't support nterms > 1 right now. I was looking at that a bit earlier, and it definitely seems like there's an opportunity for a speedup there. The Astropy implementation uses many trig_sum() calls, which finufft replaces: https://github.com/astropy/astropy/blob/bf1103e5061d4a13c928afcaa66865154f5d9afe/astropy/timeseries/periodograms/lombscargle/implementations/fastchi2_impl.py#L104. I'll open an issue to track that.

[JulienPeloton]

thanks! Let me know if you need beta testing, I would be happy to provide feedback :-)

[lgarrison]

Hi @JulienPeloton, we've implemented support for nterms > 1 in v1.1.0. It is currently a pre-release on PyPI, so you can install it with pip install --pre nifty-ls. If you would like to test it, we would appreciate any feedback you have!

[JulienPeloton]

Hi @lgarrison -- thank you for the implementation! I simply tried:

from astropy.timeseries import LombScargleMultiband
import nifty_ls  # noqa: F401

model = LombScargleMultiband(
    <time>,
    <mag>,
    <bands>,
    <errors>,
    nterms_base=2,
    nterms_band=1,
)

frequency, power = model.autopower(
    method="fast",
    sb_method="fastnifty",
    minimum_frequency=1 / period_max,
    maximum_frequency=1 / period_min,
)

But I get the following error message:

Traceback (most recent call last):
  File "/home/peloton/codes/fink-science-portal/apps/plotting.py", line 1257, in plot_variable_star
    frequency, power = model.autopower(
                       ^^^^^^^^^^^^^^^^
  File "/home/peloton/codes/fink-science-portal/portalenv/lib/python3.12/site-packages/astropy/timeseries/periodograms/lombscargle_multiband/core.py", line 381, in autopower
    power = self.power(
            ^^^^^^^^^^^
  File "/home/peloton/codes/fink-science-portal/portalenv/lib/python3.12/site-packages/astropy/timeseries/periodograms/lombscargle_multiband/core.py", line 423, in power
    power = lombscargle_multiband(
            ^^^^^^^^^^^^^^^^^^^^^^
  File "/home/peloton/codes/fink-science-portal/portalenv/lib/python3.12/site-packages/astropy/timeseries/periodograms/lombscargle_multiband/implementations/main.py", line 56, in lombscargle_multiband
    power = lombscargle_mbfast(
            ^^^^^^^^^^^^^^^^^^^
  File "/home/peloton/codes/fink-science-portal/portalenv/lib/python3.12/site-packages/astropy/timeseries/periodograms/lombscargle_multiband/implementations/mbfast_impl.py", line 29, in lombscargle_mbfast
    lombscargle(
  File "/home/peloton/codes/fink-science-portal/portalenv/lib/python3.12/site-packages/astropy/timeseries/periodograms/lombscargle/implementations/main.py", line 226, in lombscargle
    raise ValueError(
ValueError: nterms != 1 only supported with 'chi2' or 'fastchi2' methods

Apparently astropy still prevents custom sb_method to be used. Is there anything else I should add in the code?

For reference:

(portalenv) peloton@peloton:~/codes/fink-science-portal$ pip freeze | grep nifty
nifty-ls==1.1.0rc6
(portalenv) peloton@peloton:~/codes/fink-science-portal$ pip freeze | grep astropy
astropy==6.1.7

[lgarrison]

Could you try sb_method="fastnifty_chi2"?

[JulienPeloton]

thanks -- it works like a charm! I will report tomorrow on the performance in my context :-)

[lgarrison]

Great! We should probably make the sb_method="fastnifty_chi2" selection more automatic, but it's hard since we don't control the Astropy dispatch code. I'll think about it...

[JulienPeloton]

Hello @lgarrison -- I performed benchmarks on several variable star data (hundreds of measurements per star, 2 bands).

Performances

For higher order terms, using fastnifty_chi2 brings an impressive speed-up compared to the pure astropy implementation! The main driver in terms of performance is nterms_base (changing nterms_band does not change much the computation time). For illustration, here are the execution times (in seconds) for model.autopower with different backend on one typical astronomical object (~500 measurements in 2 bands):

nterms_base	nifty (s)	astropy (s)
2	0.22 +/- 0.03	6.55 +/- 0.55
3	0.31 +/- 0.03	16.63 +/- 1.02
4	0.39 +/- 0.04	25.39 +/- 1.02
5	0.47 +/- 0.06	56.11 +/- 3.29

In all cases, I am using method="fast". For astropy, using sb_method=auto or sb_method=chi2 or sb_method=fastchi2 gives the same results here. I am purely using the CPU (4 cores, 11th Gen Intel Core i7-1165G7 @ 2.80GHz), and I performed 10 loops for each cases. It is interesting to note that the speed-up increases with the number of base terms: x30, x50, x70, x120!

Accuracy

In all tested cases, the nifty and pure astropy predictions agree perfectly on finding the period. Looking closely to the residual in the periodograms (1 - nifty/astropy), it highly depends on the astropy sb_method chosen. Most users will use auto, which gives large residuals:

The residuals are several %, but tends to shrink for higher order with a slight bias (but again the best period remains the same for the two methods). If one selects the chi2 method for astropy, the residuals is much better (yscale is 1e-7):

Which means that not only nifty is faster, but also consistently accurate. Well done!