Rybicki Press algorithm: Difference between revisions

Content deleted Content added
Added figure for clarity
Changing short description from "An algorithm for inverting a matrix" to "Algorithm for inverting a matrix"
 
(30 intermediate revisions by 17 users not shown)
Line 1:
{{Short description|Algorithm for inverting a matrix}}
{{Notability|date=July 2015}}
{{Technical|date=August 2021}}
[[File:Extended_Sparse_MatrixExtended Sparse Matrix.png|thumb|Extended Sparse Matrix arising from a <math>10 \times 10</math> semi-separable matrix whose semi-separable rank is <math>4</math>.]]
The '''Rybicki–Press algorithm''' is a fast direct algorithm for inverting a matrix, whose entries are given by <math>A(i,j) = \exp(-a \vert t_i - t_j \vert)</math>, where <math>a \in \mathbb{R}</math>.<ref>{{citation
|last1 = Rybicki|first1 = George B.|last2 = Press|first2 = William H.|arxiv = comp-gas/9405004|doi = 10.1103/PhysRevLett.74.1060|journal = Physical Review Letters|page = 1060|title = Class of fast methods for processing Irregularly sampled or otherwise inhomogeneous one-dimensional data|volume = 74|year = 1995|bibcode = 1995PhRvL..74.1060R|pmid=10058924}} {{Open access}}</ref> 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.<ref>{{Cite journal|url = |title = Interpolation, realization, and reconstruction of noisy, irregularly sampled data|last = Rybicki|first = George B.|date = October 1992|journal = The Astrophysical Journal|doi = 10.1086/171845|pmid = |last2 = Press|first2 = William H.|bibcode = 1992ApJ...398..169R|volume=398|page=169}}{{Open access}}</ref><ref name=":0">{{Cite journal|last=MacLeod|first=C. L.|last2=Brooks|first2=K.|last3=Ivezic|first3=Z.|last4=Kochanek|first4=C. S.|last5=Gibson|first5=R.|last6=Meisner|first6=A.|last7=Kozlowski|first7=S.|last8=Sesar|first8=B.|last9=Becker|first9=A. C.|date=2011-02-10|title=Quasar Selection Based on Photometric Variability|url=http://arxiv.org/abs/1009.2081|journal=The Astrophysical Journal|volume=728|issue=1|pages=26|doi=10.1088/0004-637X/728/1/26|issn=0004-637X}}</ref> The most common use of the algorithm as of 2015 is in the detection of periodicity in astronomical observations<ref name=":0" />. Recently, this has been extended (Generalized Rybicki Press algorithm) for inverting matrices whose entries of the form <math>A(i,j) = \sum_{k=1}^p a_k \exp(-\beta_k \vert t_i - t_j \vert)</math><ref>{{Cite journal|last=Ambikasaran|first=Sivaram|date=2015-12-01|title=Generalized Rybicki Press algorithm|url=https://onlinelibrary.wiley.com/doi/pdf/10.1002/nla.2003|journal=Numerical Linear Algebra with Applications|language=en|volume=22|issue=6|pages=1102–1114|doi=10.1002/nla.2003|issn=1099-1506}}</ref>. The key observation for the Generalized Rybicki Press (GPP) algorithm is that the matrix <math>A</math> is a semi-separable matrix with rank <math>p</math>. More precisely, if the matrix <math>A \in \mathbb{R}^{n\times n}</math> has a semi-separable rank is <math>p</math>, the cost for solving the linear system <math>Ax=b</math> and obtain the determinant of the matrix <math>A</math>scales as <math>\mathcal{O}\left(p^2n \right)</math>, thereby making it extremely attractive for large matrices. This implementation of the GPP algorithm can be found here<ref>{{Cite web|url=https://github.com/sivaramambikasaran/ESS|title=sivaramambikasaran/ESS|website=GitHub|language=en|access-date=2018-04-05}}</ref>. The key idea is that the dense matrix <math>A</math> can be converted into a sparser matrix of a larger size (see figure on the right), whose structure can be leveraged to reduce the computational cost.
 
The '''Rybicki–Press algorithm''' is a fast [[algorithm]] for inverting a [[Matrix (mathematics)|matrix]] whose entries are given by <math>A(i,j) = \exp(-a \vert t_i - t_j \vert)</math>, where <math>a \in \mathbb{R}</math><ref name=":2">{{cite journal |last1=Rybicki |first1=George B. |last2=Press |first2=William H. |arxiv=comp-gas/9405004 |doi=10.1103/PhysRevLett.74.1060 |journal=Physical Review Letters |title=Class of fast methods for processing Irregularly sampled or otherwise inhomogeneous one-dimensional data |volume=74 |issue=7 |pages=1060–1063 |date=1995 |bibcode = 1995PhRvL..74.1060R|pmid=10058924|s2cid = 17436268}} {{Open access}}</ref> and where the <math>t_i</math> are sorted in order.<ref name=":3" /> The key observation behind the Rybicki-Press observation is that the [[matrix inverse]] of such a matrix is always a [[tridiagonal matrix]] (a matrix with nonzero entries only on the main diagonal and the two adjoining ones), and [[Tridiagonal matrix algorithm|tridiagonal systems of equations]] can be solved efficiently (to be more precise, in linear time).<ref name=":2" /> 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.<ref>{{Cite journal|title = Interpolation, realization, and reconstruction of noisy, irregularly sampled data|last1 = Rybicki|first1 = George B.|date = October 1992|journal = The Astrophysical Journal|doi = 10.1086/171845|last2 = Press|first2 = William H.|bibcode = 1992ApJ...398..169R|volume=398|page=169}}{{Open access}}</ref><ref name=":0">{{Cite journal|last1=MacLeod|first1=C. L.|last2=Brooks|first2=K.|last3=Ivezic|first3=Z.|last4=Kochanek|first4=C. S.|last5=Gibson|first5=R.|last6=Meisner|first6=A.|last7=Kozlowski|first7=S.|last8=Sesar|first8=B.|last9=Becker|first9=A. C.|date=2011-02-10|title=Quasar Selection Based on Photometric Variability|journal=The Astrophysical Journal|volume=728|issue=1|pages=26|doi=10.1088/0004-637X/728/1/26|issn=0004-637X|arxiv=1009.2081|bibcode=2011ApJ...728...26M|s2cid=28219978}}</ref> The most common use of the algorithm is in the detection of periodicity in astronomical observations{{Verify source|date=October 2021}}, such as for detecting [[Quasar|quasars]].<ref name=":0" />
The fact that matrix <math>A</math> is a semi-separable matrix also forms the basis for celerite<ref>{{Cite web|url=http://celerite.readthedocs.io/en/stable/|title=celerite — celerite 0.3.0 documentation|website=celerite.readthedocs.io|language=en|access-date=2018-04-05}}</ref> library, which is a library for fast and scalable Gaussian Process (GP) Regression in one dimension<ref name=":1">{{Cite journal|last=Foreman-Mackey|first=Daniel|last2=Agol|first2=Eric|last3=Ambikasaran|first3=Sivaram|last4=Angus|first4=Ruth|date=2017|title=Fast and Scalable Gaussian Process Modeling with Applications to Astronomical Time Series|url=http://stacks.iop.org/1538-3881/154/i=6/a=220|journal=The Astronomical Journal|language=en|volume=154|issue=6|pages=220|doi=10.3847/1538-3881/aa9332|issn=1538-3881}}</ref> with implementations in C++, Python, and Julia. The celerite method<ref name=":1" /> 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<ref>{{Cite journal|last=Foreman-Mackey|first=Daniel|date=2018|title=Scalable Backpropagation for Gaussian Processes using Celerite|url=http://stacks.iop.org/2515-5172/2/i=1/a=31|journal=Research Notes of the AAS|language=en|volume=2|issue=1|pages=31|doi=10.3847/2515-5172/aaaf6c|issn=2515-5172}}</ref><ref>{{Cite book|url=https://link.springer.com/referenceworkentry/10.1007/978-3-319-30648-3_149-1|title=Handbook of Exoplanets|last=Parviainen|first=Hannu|date=2018|publisher=Springer, Cham|isbn=9783319306483|pages=1–24|language=en|doi=10.1007/978-3-319-30648-3_149-1.pdf}}</ref>.
 
The method has been extended to the '''Generalized Rybicki-Press algorithm''' for inverting matrices with entries of the form <math>A(i,j) = \sum_{k=1}^p a_k \exp(-\beta_k \vert t_i - t_j \vert)</math>.<ref name=":3">{{Cite journal|last=Ambikasaran|first=Sivaram|date=2015-12-01|title=Generalized Rybicki Press algorithm|journal=Numerical Linear Algebra with Applications|language=en|volume=22|issue=6|pages=1102–1114|doi=10.1002/nla.2003|issn=1099-1506|arxiv=1409.7852|s2cid=1627477}}</ref> The key observation in the Generalized Rybicki-Press (GRP) algorithm is that the matrix <math>A</math> is a [[semi-separable matrix]] with rank <math>p</math> (that is, a matrix whose upper half, not including the main diagonal, is that of some matrix with [[matrix rank]] <math>p</math> and whose lower half is also that of some possibly different rank <math>p</math> matrix<ref name=":3" />) and so can be embedded into a larger [[band matrix]] (see figure on the right), whose sparsity structure can be leveraged to reduce the computational complexity. As the matrix <math>A \in \mathbb{R}^{n\times n}</math> has a semi-separable rank of <math>p</math>, the [[computational complexity]] of solving the linear system <math>Ax=b</math> or of calculating the determinant of the matrix <math>A</math> scales as <math>\mathcal{O}\left(p^2n \right)</math>, thereby making it attractive for large matrices.<ref name=":3" />
==References==
{{Reflist}}.
 
The fact that matrix <math>A</math> is a semi-separable matrix also forms the basis for {{proper name|celerite}}<ref>{{Cite web|url=httphttps://celerite.readthedocs.io/en/stable/|title=celerite — celerite 0.3.0 documentation|website=celerite.readthedocs.io|language=en|access-date=2018-04-05}}</ref> library, which is a library for fast and scalable [[Gaussian Process (GP)process Regressionregression]] in one dimension<ref name=":1">{{Cite journal|lastlast1=Foreman-Mackey|firstfirst1=Daniel|last2=Agol|first2=Eric|last3=Ambikasaran|first3=Sivaram|last4=Angus|first4=Ruth|date=2017|title=Fast and Scalable Gaussian Process Modeling with Applications to Astronomical Time Series|url=http://stacks.iop.org/1538-3881/154/i=6/a=220|journal=The Astronomical Journal|language=en|volume=154|issue=6|pages=220|doi=10.3847/1538-3881/aa9332|issn=1538-3881|arxiv=1703.09710|bibcode=2017AJ....154..220F|s2cid=88521913 |doi-access=free }}</ref> with implementations in [[C++]], [[Python (programming language)|Python]], and [[Julia (programming language)|Julia]]. The {{proper name|celerite}} method<ref name=":1" /> also provides an algorithm for generating samples from a high-dimensional distribution. The method has found attractive applications in a wide range of fields,{{Which|date=October 2021}} especially in astronomical data analysis.<ref>{{Cite journal|last=Foreman-Mackey|first=Daniel|date=2018|title=Scalable Backpropagation for Gaussian Processes using Celerite|url=http://stacks.iop.org/2515-5172/2/i=1/a=31|journal=Research Notes of the AAS|language=en|volume=2|issue=1|pages=31|doi=10.3847/2515-5172/aaaf6c|issn=2515-5172|arxiv=1801.10156|bibcode=2018RNAAS...2...31F|s2cid=102481482 |doi-access=free }}</ref><ref>{{Cite book|url=https://link.springer.com/referenceworkentry/10.1007/978-3-319-30648-3_149-1|title=Handbook of Exoplanets|last=Parviainen|first=Hannu|date=2018|publisher=Springer, Cham|isbn=9783319306483|pages=1–24|language=en|doi=10.1007/978-3-319-30648-3_149-1|chapter = Bayesian Methods for Exoplanet Science|arxiv = 1711.pdf03329}}</ref>.
[[Category:Numerical linear algebra]]
 
==See also==
* [[Invertible matrix]]
* [[Matrix decomposition]]
* [[Multidimensional signal processing]]
* [[System of linear equations]]
 
==References==
{{Signal-processing-stub}}
{{Reflist}}.
 
== External links ==
* [https://github.com/sivaramambikasaran/ESS Implementation of the Generalized Rybicki Press algorithm]
* [https://github.com/dfm/celerite celerite library on GitHub]
 
[[Category:Numerical linear algebra]]