Rybicki Press algorithm

Extended Sparse Matrix arising from a $10\times 10$ semi-separable matrix whose semi-separable rank is $4$.

The Rybicki–Press algorithm is a fast direct algorithm for inverting a matrix, whose entries are given by $A(i,j)=\exp(-a\vert t_{i}-t_{j}\vert )$, where $a\in \mathbb {R}$. It is a computational optimization of a general set of statistical methods developed to determine whether two noisy, irregularly sampled data sets are, in fact, dimensionally shifted representations of the same underlying function. The most common use of the algorithm is in the detection of periodicity in astronomical observations.

Recently, this method has been extended (Generalized Rybicki Press algorithm) for inverting matrices whose entries of the form $A(i,j)=\sum _{k=1}^{p}a_{k}\exp(-\beta _{k}\vert t_{i}-t_{j}\vert )$. The key observation in the Generalized Rybicki Press (GPP) algorithm is that the matrix $A$ is a semi-separable matrix with rank $p$. More precisely, if the matrix $A\in \mathbb {R} ^{n\times n}$ has a semi-separable rank is $p$, the cost for solving the linear system $Ax=b$ and obtaining the determinant of the matrix scales as ${\mathcal {O}}\left(p^{2}n\right)$, thereby making it extremely attractive for large matrices. This implementation of the GPP algorithm can be found here. The key idea is that the dense matrix $A$ can be converted into a sparser matrix of a larger size (see figure on the right), whose sparsity structure can be leveraged to reduce the computational complexity.

The fact that matrix $A$ is a semi-separable matrix also forms the basis for celerite library, which is a library for fast and scalable Gaussian Process (GP) Regression in one dimension with implementations in C++, Python, and Julia. The celerite method also provides an algorithm for generating samples from a high-dimensional distribution. The method has found attractive applications in a wide range of fields, especially in astronomical data analysis.

References

.