ENH: signal.lombscargle: update to the generalized Lomb-Scargle periodogram - #21277

Thanks @adammj, could you show some benchmarks to see how performance changes by removing the compiled code?

Yup. I'll get to work on some benchmarks. I only have two machines, so it'll be limited in diversity.

Benchmarks (code at the bottom)

I ran this on an M1 MacBook (macOS 14.5) and a Xeon desktop (Ubuntu 22.04), using the latest scipy (1.14.0) and numpy (2.0.1). I sampled a large range of data points (8-16384) for 3 reasonable frequency counts (512, 1024, 2048).

(The benchmark code is a little ugly and the plots could be prettier, especially the legends.)

The times, especially for a small sample sizes, are just so short that I feel the differences are basically inconsequential. Especially given the simpler code and ease of maintenance.

M1 MacBook

Xeon desktop

import numpy as np
import seaborn as sns
from scipy.signal import lombscargle
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import figure_standard  # a library of mine for consistent plots


figure_standard.standardize_plot_settings()
rng = np.random.default_rng(12)
A = 2.1
w0 = 1.12  # rad/sec

gs = gridspec.GridSpec(1, 3)
ax = []
for i in range(3):
    ax += [plt.subplot(gs[0, i])]



for i, n_freqs_i in enumerate([9, 10, 11]):
    n_freqs = 2**n_freqs_i
    print(f"n_freqs: {n_freqs}")

    time_n = []
    old_avg = []
    old_std = []
    new_avg = []
    new_std = []

    for j in range(3, 15):
        n_times = 2**j
        time_n += [n_times]
        
        t = rng.uniform(0, 10*np.pi, n_times)
        y = A * np.cos(w0*t)
        w = np.linspace(0.01, 10, n_freqs)

        def test_old():
            pgram = lombscargle(t, y, w, normalize=True)  # existing scipy function

        def test_new():
            pgram_new = lombscargle_new(t, y, w, normalize=True)  # just the new function

        result_old = %timeit -n 9 -r 9 -o test_old()
        old_avg += [result_old.average]
        old_std += [result_old.stdev]

        result_new = %timeit -n 9 -r 9 -o test_new()
        new_avg += [result_new.average]
        new_std += [result_new.stdev]

    time_np = np.array(time_n)
    old_avg_np = np.array(old_avg)/9
    new_avg_np = np.array(new_avg)/9
    old_std_np = np.array(old_std)/9
    new_std_np = np.array(new_std)/9

    sns.lineplot(x=time_np, y=old_avg_np, ax=ax[i], legend=True)
    ax[i].fill_between(time_np, old_avg_np-old_std_np, old_avg_np+old_std_np, alpha=0.2)

    sns.lineplot(x=time_np, y=new_avg_np, ax=ax[i], legend=True)
    ax[i].fill_between(time_np, new_avg_np-new_std_np, new_avg_np+new_std_np, alpha=0.2)

    ax[i].set_ylim([0.00001, 1])
    ax[i].set_xlim([5, 20000])
    ax[i].set_yscale('log')
    ax[i].set_xscale('log')

    if i == 0:
        ax[i].legend(labels=['current',None,None,'new'], loc='upper left')
        ax[i].set_ylabel('time / iteration (sec)')
    
    ax[i].set_xlabel(f'# data points (for {n_freqs} freqs)')


plt.tight_layout()
plt.subplots_adjust(wspace=0.25, hspace=0)
figure_standard.scale_figure_by_axes(ax[0], 1.5, 1.5)


plt.savefig("result_linux.png",
            bbox_inches="tight",
            pad_inches=0.05,
            transparent=False,
        )

Thanks for the benchmarks, I am not a signal person so don't have a sense of whether that performance hit is acceptable so I will leave that decision to another maintainer.

I finally found some time to familiarize myself with the (generalized) Lomb–Scargle periodogram a bit. I think this PR would be a valuable addition—especially, being able to also determine the spectrum's phase is a big plus, IMO. I do have a couple of remarks:

• The docstring could use some love:

  • It would be nice, if the documentation stated that this function is based on a weighted least-squares fit of $y_i = a\cos(\omega t_i) + b\cos(\omega t_i) + o_i$ (something similar like in Wikipedia). An advantage to a standard least squares calculation, as in scipy.linalg.lstsq(), is that this implementation is typically numerically more stable.
  • It is important to mention that the spectral values are calculated independently of each other (due to the for loop in the code). This is a difference to the more common FFT based methods. Especially in the presence of noise, the result differ quite a bit to case were all parameters would be calculated at once (e.g., there may a different offset for each frequency).
  • The Lomb-Scargle Periodograms section in the SciPy User Guide is probably not useful anymore and could be removed.

• I am not worried much of the slowdown of the new version, since for more data points, the speed difference diminishes according to the nice benchmarks.

• As far as I can tell, the unit tests do not check, if the phase angle is calculated correctly for normalize == 'amplitude'.

• Keeping the classic version, i.e., not determining the offset (precenter == False) is still useful for some fringe applications: Imagine having data points from a sine wave with known frequency (no offset) and you are interested in amplitude and phase. If the data points only span a (small) part of the cycle period, subtracting the mean would be counterproductive.

@DietBru Thank you! I'll get to work on your suggestions.

On this point, however, what are you envisioning?

• Keeping the classic version, i.e., not determining the offset (precenter == False) is still useful for some fringe applications: Imagine having data points from a sine wave with known frequency (no offset) and you are interested in amplitude and phase. If the data points only span a (small) part of the cycle period, subtracting the mean would be counterproductive.

I'll have to give this some more thought, but I can think of a few issues:

1. Regardless of what the lowest frequency is, if the total time is only equal to half a cycle of that lowest frequency, then both formulations produce not very useful results for that frequency. The peak becomes -very- broad (as in, it's looks more like a rolling hill, than a spectral peak). I can demonstrate this with code+figure.
2. The generalized version is supposed to catch any unknown offset regardless of if the user knows it exists. To me, trying to prevent this behavior would at the very least strongly suggest that the parameter default should be also changed to precenter=True. Because the default behavior, I believe, should be to make use of the generalized version's improvements.
3. Alternatively, add a new parameter. The problem is the name precenter gives the impression that you're performing the operation y = y - y.mean() first. When the generalized version is not doing that at all. So, if the user explicitly wants the code to assume offset=0 for all frequencies, I think that an additional new parameter name like floating_mean (with a default =True) would more appropriately match the intention of what you're after, and the standard nomenclature of what the generalized version is doing.
4. I'll have to think through how best I could "add back" a no-offset version. I'm thinking it would require a conditional outside the for loop (as in, two separate for loops, one for each branch), because all of the equations inside the loop are predicated on the generalized version (and adding a conditional inside the loop would be a huge performance penalty).

Update

Well, that was easier than I thought. I added a new parameter, floating_mean (default = True). Only one conditional is required inside the for loop, and another inside the calculation for the normalize=='normalize' condition. With only a single additional conditional per iteration, there won't be any performance issue.

# inside for loop
if floating_mean:
    # calculate best-fit offset for each frequency independently (default)
    C_sum = wcoswt.sum()
    S_sum = wsinwt.sum()
    YC = YC_hat - Y_sum * C_sum
    YS = YS_hat - Y_sum * S_sum
    CC = CC_hat - C_sum * C_sum
    SS = SS_hat - S_sum * S_sum
    CS = CS_hat - C_sum * S_sum
else:
    # assume offset=0 for all frequencies
    YC = YC_hat
    YS = YS_hat
    CC = CC_hat
    SS = SS_hat
    CS = CS_hat

D = CC * SS - CS * CS

# inside if normalize == "normalize" condition
if floating_mean:
    YY: float = YY_hat - Y_sum * Y_sum
else:
    YY: float = YY_hat

@DietBru

You said:

• The Lomb-Scargle Periodograms section in the SciPy User Guide is probably not useful anymore and could be removed.

I think it should stick around. Over the years, when discussing lombscargle, I've had to spend immense amount of time explaining LSSA, and how, mathematically it reduces to the same thing as Fourier analysis when using uniformly sampled data. (Until then, they think it's some "fake" frequency analysis.)

However, I think everything after and including "The lombscargle function calculates the periodogram using a slightly modified algorithm due to Townsend..." should be updated to mention the generalized version. Furthermore, the equations at the end should be removed, as I don't think think they add anything to the discussion.

FYI I have identified this patch as the root cause of some new failures in astropy's test suite: astropy/astropy#17149

I haven't scanned through the whole discussion but since the OP opens as "The function is a backwards-compatible replacement" I'd like to ask first: is it expected to produce (even slightly) different results ? I'll have more time to go through the whole history tomorrow, apologies if my question is already answered somewhere in the discussion.

Thanks for explaining, this really isn't my domain of expertise so do you mind clarifying what you mean by better? Does this mean more accurate?

For data that has an unknown offset and especially data that has clipped peaks or troughs, yes, it is more accurate. For "simple" data, the results will be nearly identical (except, again, for the numerical differences inherent in floating point math). But there are additional features (like weights, and being able to produce the complex amplitude) that the old formulation did not allow for.

My question about backwards compatibility was mainly about the stronger input validation, are there any inputs that previously produced a valid output that might now raise an error?

Thanks for reporting @neutrinoceros this should be the relevant part of the discussion. However please do let us know if you find the above not to be the case.

I'd like to ask first: is it expected to produce (even slightly) different results ?

@neutrinoceros Yes, the values posted in your downstream issue in astropy/astropy#17149 just highlight numerical differences that are always present when using floating-point math. Even if two formulas mathematically reduce to the same exact thing, the path your calculations take using float32, float64, floatWHATEVER, will wind up with ever so slightly different results.

Thanks for your prompt responses !

The differences we get seem kinda high for numerical noise, for instance

Max absolute difference among violations: 0.05446034
Max relative difference among violations: 4.22703988

nonetheless, those are not necessarily direct outputs from scipy. Further inquiry is needed on the side of astropy for now.

It looks like the issue was that it is not really fully backwards-compatible with the default of floating_mean=True, setting this to False seems to completely restore our previous results. Still need to discuss if the new defaults are not actually preferable, though.

@dhomeier Defaulting to using the floating mean formulation (per Zechmeister and Kürster) was discussed above, as the previous precenter parameter (which was kept for backwards compatibility of the API) is more misleading than correct--especially in astronomy light curves, where observation times might systematically prevent observing the peaks and/or troughs, which will bias the mean.

@adammj yes, I was assuming something like that, just have not had the time to read up on the discussion.
I assumed our own pure Python + Numpy implementation must have had some motivation to use the previous convention, and in fact it does refer to Zechmeister & Kürster, though "using the floating mean" seems to be understood in the opposite sense there.
I shall propose that we at least allow using floating_mean=True in new releases, and if possible add a corresponding option to our internal implementation.

I shall propose that we at least allow using floating_mean=True in new releases, and if possible add a corresponding option to our internal implementation.

Ah, nevermind – looks like this is exactly what @jakevdp has already implemented as fit_mean on our side (and was disabled for the scipy method ;-), so probably all that needs to be done is adding this to our scipy wrapper when available.

The PR description notes backward compatibility, but it should be updated to make clear that although the API is backward compatible, the outputs are not: the function outputs a different statistic than it did previously. Namely, fit_mean is set to True by default, and to recover previous outputs it must explicitly be set to False.

I also think it would be worth discussing whether this breakage is advisable. I'd lean toward fit_mean=False by default (at least for now) to avoid this kind of breakage, despite the fact that all else being equal, fit_mean=True would be a better default.

@tylerjereddy I wasn't sure who to notify, but I see your username next to each of the recent releases. As @mhvk mentions in astropy/astropy#17149, it's worth highlighting this change in the release notes. While the API is backwards compatible, it may cause some results to change with respect to the previous version's output, if the data has a vertical offset (even if it is removed with the precenter=True argument). Any differences will be "for the better", as in, the new implementation is more appropriately handling these offsets by default. But, if they are comparing results between versions, and had an offset, they may see a difference.

@tylerjereddy I wasn't sure who to notify, but I see your username next to each of the recent releases. As @mhvk mentions in astropy/astropy#17149, it's worth highlighting this change in the release notes. While the API is backwards compatible, it may cause some results to change with respect to the previous version's output, if the data has a vertical offset (even if it is removed with the precenter=True argument). Any differences will be "for the better", as in, the new implementation is more appropriately handling these offsets by default. But, if they are comparing results between versions, and had an offset, they may see a difference.

Thanks for explaining. I think I'd prefer to be on the conservative side, as unexpected changes to values can cause headaches for our downstream users. Perhaps we should set fit_mean=False with a future warning that it will change to true?

Thanks for explaining. I think I'd prefer to be on the conservative side, as unexpected changes to values can cause headaches for our downstream users. Perhaps we should set fit_mean=False with a future warning that it will change to true?

@j-bowhay I think the previous 3 comments by @jakevdp and @dhomeier were supposed to go on astropy/astropy#17149, as they're talking about astropy features. Specifically, fit_mean is an internal astropy implementation detail, and not part of scipy.signal.lombscargle.

Thanks for explaining. I think I'd prefer to be on the conservative side, as unexpected changes to values can cause headaches for our downstream users. Perhaps we should set fit_mean=False with a future warning that it will change to true?

@j-bowhay I think the previous 3 comments by @jakevdp and @dhomeier were supposed to go on astropy/astropy#17149, as they're talking about astropy features. Specifically, fit_mean is an internal astropy implementation detail, and not part of scipy.signal.lombscargle.

Ah thanks, got confused trying to follow on my phone

@adammj I suggested keeping the more conservative argument in the cross-linked issue and issuing a warning. If the astropy folks are actually happy in the end though it may be "ok," I'll ping them over there.

One place where fit_mean=True will cause problems is when you use lombscargle to find the window frequency by passing in a unit signal, like this:

signal_psd = lombscargle(x, y, freqs)
window_psd = lombscargle(x, np.ones_like(y), freqs)

If fit_mean is True by default, then window_psd is under-constrained, and you'd have to explicitly set fit_mean=False to get meaningful output. This is not "backward-compatible" even in the looser sense of the term.

@jakevdp I'm confused. fit_mean is a astropy.lombscargle parameter, and it is true by default. And, as best as I can tell, it is equivalent (in spirit) to floating_mean in scipy.lombscargle, which (currently) is also true by default. So, it looks like to me, the two methods are in agreement that a "floating mean" is best, and should be the default.

Yes, in astropy it has always been set to True, because it's the more sensible default when the option is available.

Here the argument is that changing the default is a breaking change, and so it shouldn't be done quietly because it will break downstream projects (as it did in the linked astropy issue).

Moving discussion to #21696 (comment)