The Lomb-Scargle Periodogram is a well-known method of finding periodicity in irregularly-sampled time-series data.
The common implementation of the periodogram is relatively slow: for $N$ data points, a frequency grid of $\sim N$ frequencies is required and the computation scales as $O[N^2]$.
In a 1989 paper, Press and Rybicki presented a faster technique which makes use of fast Fourier transforms to reduce this cost to $O[N\log N]$ on a regular frequency grid.
The `gatspy`

package implement this in the `LombScargleFast`

object, which we'll explore below.

But first, we'll motivate *why* this algorithm is needed at all.
We'll start this notebook with some standard imports:

In [1]:

```
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
# use seaborn's default plotting styles for matplotlib
import seaborn; seaborn.set()
```

In [2]:

```
def create_data(N, period=2.5, err=0.1, rseed=0):
rng = np.random.RandomState(rseed)
t = np.arange(N, dtype=float) + 0.3 * rng.randn(N)
y = np.sin(2 * np.pi * t / period) + err * rng.randn(N)
return t, y, err
t, y, dy = create_data(100, period=20)
plt.errorbar(t, y, dy, fmt='o');
```

From this, our algorithm should be able to identify any periodicity that is present.

The Lomb-Scargle Periodogram works by evaluating a power for a set of candidate frequencies $f$. So the first question is, how many candidate frequencies should we choose?
It turns out that this question is *very* important. If you choose the frequency spacing poorly, it may lead you to miss strong periodic signal in the data!

First, let's think about the frequency spacing we need in our grid. If you're asking about a candidate frequency $f$, then data with range $T$ contains $T \cdot f$ complete cycles. If our error in frequency is $\delta f$, then $T\cdot\delta f$ is the error in number of cycles between the endpoints of the data. If this error is a significant fraction of a cycle, this will cause problems. This givs us the criterion $$ T\cdot\delta f \ll 1 $$ Commonly, we'll choose some oversampling factor around 5 and use $\delta f = (5T)^{-1}$ as our frequency grid spacing.

Next, we need to choose the limits of the frequency grid. On the low end, $f=0$ is suitable, but causes some problems – we'll go one step away and use $\delta f$ as our minimum frequency. But on the high end, we need to make a choice: what's the highest frequency we'd trust our data to be sensitive to? At this point, many people are tempted to mis-apply the Nyquist-Shannon sampling theorem, and choose some version of the Nyquist limit for the data. But this is entirely wrong! The Nyquist frequency applies for regularly-sampled data, but irregularly-sampled data can be sensitive to much, much higher frequencies, and the upper limit should be determined based on what kind of signals you are looking for.

Still, a common (if dubious) rule-of-thumb is that the high frequency is some multiple of what Press & Rybicki call the "average" Nyquist frequency, $$ \hat{f}_{Ny} = \frac{N}{2T} $$ With this in mind, we'll use the following function to determine a suitable frequency grid:

In [3]:

```
def freq_grid(t, oversampling=5, nyquist_factor=3):
T = t.max() - t.min()
N = len(t)
df = 1. / (oversampling * T)
fmax = 0.5 * nyquist_factor * N / T
N = int(fmax // df)
return df + df * np.arange(N)
```

Now let's use the `gatspy`

tools to plot the periodogram:

In [4]:

```
t, y, dy = create_data(100, period=2.5)
freq = freq_grid(t)
print(len(freq))
```

750

In [5]:

```
from gatspy.periodic import LombScargle
model = LombScargle().fit(t, y, dy)
period = 1. / freq
power = model.periodogram(period)
plt.plot(period, power)
plt.xlim(0, 5);
```

The algorithm finds a strong signal at a period of 2.5.

To demonstrate explicitly that the Nyquist rate doesn't apply in irregularly-sampled data, let's use a period below the averaged sampling rate and show that we can find it:

In [6]:

```
t, y, dy = create_data(100, period=0.3)
period = 1. / freq_grid(t, nyquist_factor=10)
model = LombScargle().fit(t, y, dy)
power = model.periodogram(period)
plt.plot(period, power)
plt.xlim(0, 1);
```

With these rules in mind, we see that the size of the frequency grid is approximately $$ N_f = \frac{f_{max}}{\delta f} \propto \frac{N/(2T)}{1/T} \propto N $$ So for $N$ data points, we will require some multiple of $N$ frequencies (with a constant of proportionality typically on order 10) to suitably explore the frequency space. This is the source of the $N^2$ scaling of the typical periodogram: finding periods in $N$ datapoints requires a grid of $\sim 10N$ frequencies, and $O[N^2]$ operations. When $N$ gets very, very large, this becomes a problem.

`LombScargleFast`

¶Finally we get to the meat of this discussion.

In a 1989 paper, Press and Rybicki proposed a clever method whereby a Fast Fourier Transform is used on a grid *extirpolated* from the original data, such that this problem can be solved in $O[N\log N]$ time. The `gatspy`

package contains a pure-Python implementation of this algorithm, and we'll explore it here.
If you're interested in seeing how the algorithm works in Python, check out the code in the gatspy source.
It's far more readible and understandable than the Fortran source presented in Press *et al.*

For convenience, the implementation has a `periodogram_auto`

method which automatically selects a frequency/period range based on an oversampling factor and a nyquist factor:

In [7]:

```
from gatspy.periodic import LombScargleFast
help(LombScargleFast.periodogram_auto)
```

In [8]:

```
from gatspy.periodic import LombScargleFast
t, y, dy = create_data(100)
model = LombScargleFast().fit(t, y, dy)
period, power = model.periodogram_auto()
plt.plot(period, power)
plt.xlim(0, 5);
```

`LombScargleAstroML`

(a fast implementation of the $O[N^2]$ algorithm) and `LombScargleFast`

, which is the fast FFT-based implementation:

In [9]:

```
from time import time
from gatspy.periodic import LombScargleAstroML, LombScargleFast
def get_time(N, Model):
t, y, dy = create_data(N)
model = Model().fit(t, y, dy)
t0 = time()
model.periodogram_auto()
t1 = time()
result = t1 - t0
# for fast operations, we should do several and take the median
if result < 0.1:
N = min(50, 0.5 / result)
times = []
for i in range(5):
t0 = time()
model.periodogram_auto()
t1 = time()
times.append(t1 - t0)
result = np.median(times)
return result
N_obs = list(map(int, 10 ** np.linspace(1, 4, 5)))
times1 = [get_time(N, LombScargleAstroML) for N in N_obs]
times2 = [get_time(N, LombScargleFast) for N in N_obs]
```

In [10]:

```
plt.loglog(N_obs, times1, label='Naive Implmentation')
plt.loglog(N_obs, times2, label='FFT Implementation')
plt.xlabel('N observations')
plt.ylabel('t (sec)')
plt.legend(loc='upper left');
```