{{ message }}
Repository navigation
LombScargle methods disagree, and chi2 raises LinAlgError, at frequencies where the sine term is zero #20571
Description
Description
While testing an experimental version of Perch, we found the following and figured it was worth a report.
When every sample time puts the sine term at zero for a frequency, the exact LombScargle methods stop agreeing. For example, that happens for evenly spaced integer times at frequency 1. At those frequencies:
chi2returns 0,slowreturnsinf, andcythonreturnsnanfor the same input.- At frequency 0.5, the three methods return three different values.
- With half-step times,
chi2raisesnumpy.linalg.LinAlgError: Singular matrixinstead of returning a value. It solves its least-squares problem withnp.linalg.solve, and the matrix has a zero row there.
At other frequencies all three methods agree.
Expected behavior
The methods return the same power for the same input, or all raise the same clear error. chi2 shouldn't fail with a raw LinAlgError.
How to Reproduce
- Install the development version of astropy.
- Run:
import numpy as np
from astropy.timeseries import LombScargle
t = np.arange(5.0)
y = np.array([1.0, 2.0, 1.5, 2.5, 1.0])
frequency = np.array([0.2, 0.5, 1.0])
for method in ("chi2", "slow", "cython"):
print(method, LombScargle(t, y).power(frequency, method=method))
# chi2 [0.43423329 0.84313725 0. ]
# slow [0.43423329 0.97521541 inf]
# cython [0.43423329 0.36000577 nan]
t = np.arange(0, 10, 0.5)
y = np.sin(t)
LombScargle(t, y).power(np.array([2.0]), method="chi2")
# numpy.linalg.LinAlgError: Singular matrixOn astropy 8.0.1, cython gives 1.179 at frequency 0.5, which is above 1, the upper limit for the standard normalization.
Possibly related: #11302, where fast disagreed with the other methods.
Versions
platform.platform() = 'macOS-26.2-arm64-arm-64bit'
platform.python_version() = '3.12.12'
astropy 8.1.0.dev724+g3265cc615
numpy 2.5.3
pyerfa 2.0.1.5
Activity
pllim commentedon Oct 9, 2026
I see that you mentioned #11302 , which was supposedly addressed by #17842 . Hmm, @kkrutkowski and @astrofrog ?
kkrutkowski commentedon Oct 9, 2026
I'll look into that with more detail in nearby future, but that's not related to the issue #11302 (caused by instability of well-conditioned systems) and the
# numpy.linalg.LinAlgError: Singular matrixerror as well as the nan output are correct ones I guess. The periodogram itself relies on the system of equations itself being overdetermined - that one is not.
Essentially it's the issue #10070 opened once again (exactly the same divison by zero issue). Here the result results from undefined behavior happening within C for the Cython code, although I have no idea why the
slowimplementation returnsinfinstead of the expectednan. It must be some round-off error causing non-zero numerator, but I have no idea where exactly it happens.Now - we could theoretically fix this by implementing a condition number threshold and there are a few ways to do that, however estimation of the condition number can be pretty computationally expensive and/or inaccurate.
A few methods I've tested but are mostly implementable only in C for multiharmonic case, rather, than Python can be found in https://github.com/kkrutkowski/Toeplitz-LS/blob/main/src/linalg.c. In case of the single-harmonic case (
slow,fastandcython) the codition number reduces to a pretty straightforward formula, so we would have to agree how we report singular frequencies.Reporting an error for an entire grid is generally a bad idea,
nan's can crash some NumPy/SciPy functions and generally every solution has its own downsides.Nifty-LS does considers returning using masked arrays instead of standard ones for the same reason, and I guess that could be done for AstroPy as well, but I have no idea how would it affect the backwards compatibility.
In general it's feasible to fix, but that bug is known for over six years already and I think if it was applied it would be most reasonable to do this with major solver updates (as the linear equations can be solved by applying Szego / Levinson recursion as well, with no issues) along with simple condition number fixes. Alternatively it could be applied only for single harmonic while leaving
chi²for later, which would be a pretty straightforward fix I guess.pllim commentedon Oct 9, 2026
Thanks for the clarification! So should I instead close this as a duplicate of #10070 ?
kkrutkowski commentedon Oct 9, 2026
Honestly I have no idea. Both are a different symptoms of the same undefined behavior, which is already known for years, but the symptoms themselves are different enough it might be worth leaving open, until it gets resolved in the future.
kkrutkowski commentedon Oct 10, 2026
I've just checkec the formulae. The 2-norm condition number is given by
for
fit_mean=false, andfor
fit_mean=true. In the given example values of κ2 exceed 10^30, which lies significantly above the inverse of the machine epsilon at double precision, essentially stating, that the output mush be garbage in this exact case.Now the question which remains is what we do with these. If we want to minimize potential user impact, we could just set the rejected frequencies to 0, which should be pretty straightforward. However the question about rejection criterion remains open. We could base it on:
Additionally the
chi2methods remain slightly different issue here.We could modify the standard
chi2method to use Szego recursion (more computationally efficient and numerically stable, than the current one) and then get the pretty conservative condition number estimate from the formula established by Cybenko, which is pretty conservative but usually works really well.In this case we could also easily modify the Cython implementation to accept optional
ntermsparameter as an input, since it's already written differently from the rest of the implementations and adding support for higher number of terms could be added pretty much seamlessly. Also it could potentially outperform currentfastchi2, in case of which the NumPy's solver is a pretty significant bottleneck, already for 2-term periodogram.And finally there the
fastchi2for which there is no easy way out. We could usenumpy.linalg.condif we do not want to make any significant changes, but this one is really slow. Alternatively we could rewrite a singificant part of the code in C, which could make it faster, but probably would require a lot of effort.astrofrog commentedon Oct 10, 2026
We should keep this one open (and both ideally) to not lose all the useful discussion here