HOME

TheInfoList



OR:

Levinson recursion or Levinson–Durbin recursion is a procedure in
linear algebra Linear algebra is the branch of mathematics concerning linear equations such as: :a_1x_1+\cdots +a_nx_n=b, linear maps such as: :(x_1, \ldots, x_n) \mapsto a_1x_1+\cdots +a_nx_n, and their representations in vector spaces and through matrice ...
to recursively calculate the solution to an equation involving a
Toeplitz matrix In linear algebra, a Toeplitz matrix or diagonal-constant matrix, named after Otto Toeplitz, is a matrix in which each descending diagonal from left to right is constant. For instance, the following matrix is a Toeplitz matrix: :\qquad\begin a & b ...
. The
algorithm In mathematics and computer science, an algorithm () is a finite sequence of rigorous instructions, typically used to solve a class of specific problems or to perform a computation. Algorithms are used as specifications for performing ...
runs in time, which is a strong improvement over Gauss–Jordan elimination, which runs in Θ(''n''3). The Levinson–Durbin algorithm was proposed first by Norman Levinson in 1947, improved by James Durbin in 1960, and subsequently improved to and then multiplications by W. F. Trench and S. Zohar, respectively. Other methods to process data include
Schur decomposition In the mathematical discipline of linear algebra, the Schur decomposition or Schur triangulation, named after Issai Schur, is a matrix decomposition. It allows one to write an arbitrary complex square matrix as unitarily equivalent to an upper t ...
and Cholesky decomposition. In comparison to these, Levinson recursion (particularly split Levinson recursion) tends to be faster computationally, but more sensitive to computational inaccuracies like
round-off error A roundoff error, also called rounding error, is the difference between the result produced by a given algorithm using exact arithmetic and the result produced by the same algorithm using finite-precision, rounded arithmetic. Rounding errors are d ...
s. The Bareiss algorithm for Toeplitz matrices (not to be confused with the general Bareiss algorithm) runs about as fast as Levinson recursion, but it uses space, whereas Levinson recursion uses only ''O''(''n'') space. The Bareiss algorithm, though, is numerically stable, whereas Levinson recursion is at best only weakly stable (i.e. it exhibits numerical stability for
well-conditioned In numerical analysis, the condition number of a function measures how much the output value of the function can change for a small change in the input argument. This is used to measure how sensitive a function is to changes or errors in the input ...
linear systems). Newer algorithms, called ''asymptotically fast'' or sometimes ''superfast'' Toeplitz algorithms, can solve in for various ''p'' (e.g. ''p'' = 2, ''p'' = 3 ). Levinson recursion remains popular for several reasons; for one, it is relatively easy to understand in comparison; for another, it can be faster than a superfast algorithm for small ''n'' (usually ''n'' < 256).


Derivation


Background

Matrix equations follow the form : \mathbf M \, \vec x = \vec y. The Levinson–Durbin algorithm may be used for any such equation, as long as M is a known
Toeplitz matrix In linear algebra, a Toeplitz matrix or diagonal-constant matrix, named after Otto Toeplitz, is a matrix in which each descending diagonal from left to right is constant. For instance, the following matrix is a Toeplitz matrix: :\qquad\begin a & b ...
with a nonzero main diagonal. Here \vec y is a known vector, and \vec x is an unknown vector of numbers ''x''''i'' yet to be determined. For the sake of this article, ''ê''''i'' is a vector made up entirely of zeroes, except for its ''i''th place, which holds the value one. Its length will be implicitly determined by the surrounding context. The term ''N'' refers to the width of the matrix above – M is an ''N''×''N'' matrix. Finally, in this article, superscripts refer to an ''inductive index'', whereas subscripts denote indices. For example (and definition), in this article, the matrix T''n'' is an ''n''×''n'' matrix that copies the upper left ''n''×''n'' block from M – that is, ''T''''n''''ij'' = ''M''''ij''. T''n'' is also a Toeplitz matrix, meaning that it can be written as : \mathbf T^n = \begin t_0 & t_ & t_ & \dots & t_ \\ t_1 & t_0 & t_ & \dots & t_ \\ t_2 & t_1 & t_0 & \dots & t_ \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ t_& t_ & t_ & \dots & t_0 \end.


Introductory steps

The algorithm proceeds in two steps. In the first step, two sets of vectors, called the ''forward'' and ''backward'' vectors, are established. The forward vectors are used to help get the set of backward vectors; then they can be immediately discarded. The backwards vectors are necessary for the second step, where they are used to build the solution desired. Levinson–Durbin recursion defines the ''n''th "forward vector", denoted \vec f^n, as the vector of length ''n'' which satisfies: :\mathbf T^n \vec f^n = \hat e_1. The ''n''th "backward vector" \vec b^n is defined similarly; it is the vector of length ''n'' which satisfies: :\mathbf T^n \vec b^n = \hat e_n. An important simplification can occur when ''M'' is a
symmetric matrix In linear algebra, a symmetric matrix is a square matrix that is equal to its transpose. Formally, Because equal matrices have equal dimensions, only square matrices can be symmetric. The entries of a symmetric matrix are symmetric with ...
; then the two vectors are related by ''b''''n''''i'' = ''f''''n''''n''+1−''i''—that is, they are row-reversals of each other. This can save some extra computation in that special case.


Obtaining the backward vectors

Even if the matrix is not symmetric, then the ''n''th forward and backward vector may be found from the vectors of length ''n'' − 1 as follows. First, the forward vector may be extended with a zero to obtain: :\mathbf T^n \begin \vec f^ \\ 0 \\ \end = \begin \ & \ & \ & t_ \\ \ & \mathbf T^ & \ & t_ \\ \ & \ & \ & \vdots \\ t_ & t_ & \dots & t_0 \\ \end \begin \ \\ \vec f^ \\ \ \\ 0 \\ \ \\ \end = \begin 1 \\ 0 \\ \vdots \\ 0 \\ \varepsilon_f^n \end. In going from ''T''''n''−1 to ''Tn'', the extra ''column'' added to the matrix does not perturb the solution when a zero is used to extend the forward vector. However, the extra ''row'' added to the matrix ''has'' perturbed the solution; and it has created an unwanted error term ''ε''''f'' which occurs in the last place. The above equation gives it the value of: : \varepsilon_f^n \ = \ \sum_^ \ M_ \ f_^ \ = \ \sum_^ \ t_ \ f_^. This error will be returned to shortly and eliminated from the new forward vector; but first, the backwards vector must be extended in a similar (albeit reversed) fashion. For the backwards vector, : \mathbf T^n \begin 0 \\ \vec b^ \\ \end = \begin t_0 & \dots & t_ & t_ \\ \vdots & \ & \ & \ \\ t_ & \ & \mathbf T^ & \ \\ t_ & \ & \ & \end \begin \ \\ 0 \\ \ \\ \vec b^ \\ \ \\ \end = \begin \varepsilon_b^n \\ 0 \\ \vdots \\ 0 \\ 1 \end. As before, the extra column added to the matrix does not perturb this new backwards vector; but the extra row does. Here we have another unwanted error ''ε''''b'' with value: : \varepsilon_b^n \ = \ \sum_^n \ M_ \ b_^ \ = \ \sum_^ \ t_ \ b_i^. \ These two error terms can be used to form higher-order forward and backward vectors described as follows. Using the linearity of matrices, the following identity holds for all (\alpha,\beta): :\mathbf T \left( \alpha \begin \vec f \\ \ \\ 0 \\ \end + \beta \begin 0 \\ \ \\ \vec b \end \right ) = \alpha \begin 1 \\ 0 \\ \vdots \\ 0 \\ \varepsilon_f \\ \end + \beta \begin \varepsilon_b \\ 0 \\ \vdots \\ 0 \\ 1 \end. If ''α'' and ''β'' are chosen so that the right hand side yields ''ê''1 or ''ê''''n'', then the quantity in the parentheses will fulfill the definition of the ''n''th forward or backward vector, respectively. With those alpha and beta chosen, the vector sum in the parentheses is simple and yields the desired result. To find these coefficients, \alpha^n_, \beta^n_ are such that : : \vec f^n = \alpha^n_ \begin \vec f^\\ 0 \end +\beta^n_\begin0\\ \vec b^ \end and respectively \alpha^n_, \beta^n_ are such that : :\vec b^n = \alpha^n_ \begin \vec f^\\ 0 \end +\beta^n_\begin 0\\ \vec b^ \end. By multiplying both previous equations by ^n one gets the following equation: : \begin 1 & \varepsilon^n_b \\ 0 & 0 \\ \vdots & \vdots \\ 0 & 0 \\ \varepsilon^n_f & 1 \end \begin \alpha^n_f & \alpha^n_b \\ \beta^n_f & \beta^n_b \end = \begin 1 & 0 \\ 0 & 0 \\ \vdots & \vdots \\ 0 & 0 \\ 0 & 1 \end. Now, all the zeroes in the middle of the two vectors above being disregarded and collapsed, only the following equation is left: : \begin 1 & \varepsilon^n_b \\ \varepsilon^n_f & 1 \end \begin \alpha^n_f & \alpha^n_b \\ \beta^n_f & \beta^n_b \end = \begin 1 & 0 \\ 0 & 1 \end. With these solved for (by using the Cramer 2×2 matrix inverse formula), the new forward and backward vectors are: : \vec f^n = \begin \vec f^ \\ 0 \end - \begin 0 \\ \vec b^ \end : \vec b^n = \begin 0 \\ \vec b^ \end - \begin \vec f^ \\ 0 \end. Performing these vector summations, then, gives the ''n''th forward and backward vectors from the prior ones. All that remains is to find the first of these vectors, and then some quick sums and multiplications give the remaining ones. The first forward and backward vectors are simply: : \vec f^1 = \vec b^1 = \left \right= \left \right


Using the backward vectors

The above steps give the ''N'' backward vectors for ''M''. From there, a more arbitrary equation is: : \vec y = \mathbf M \ \vec x. The solution can be built in the same recursive way that the backwards vectors were built. Accordingly, \vec x must be generalized to a sequence of intermediates \vec x^n, such that \vec x^N = \vec x. The solution is then built recursively by noticing that if : \mathbf T^ \begin x_1^ \\ x_2^ \\ \vdots \\ x_^ \\ \end = \begin y_1 \\ y_2 \\ \vdots \\ y_ \end. Then, extending with a zero again, and defining an error constant where necessary: : \mathbf T^ \begin x_1^ \\ x_2^ \\ \vdots \\ x_^ \\ 0 \end = \begin y_1 \\ y_2 \\ \vdots \\ y_ \\ \varepsilon_x^ \end. We can then use the ''n''th backward vector to eliminate the error term and replace it with the desired formula as follows: : \mathbf T^ \left ( \begin x_1^ \\ x_2^ \\ \vdots \\ x_^ \\ 0 \\ \end + (y_n - \varepsilon_x^) \ \vec b^n \right ) = \begin y_1 \\ y_2 \\ \vdots \\ y_ \\ y_n \end. Extending this method until ''n'' = ''N'' yields the solution \vec x. In practice, these steps are often done concurrently with the rest of the procedure, but they form a coherent unit and deserve to be treated as their own step.


Block Levinson algorithm

If ''M'' is not strictly Toeplitz, but block Toeplitz, the Levinson recursion can be derived in much the same way by regarding the block Toeplitz matrix as a Toeplitz matrix with matrix elements (Musicus 1988). Block Toeplitz matrices arise naturally in signal processing algorithms when dealing with multiple signal streams (e.g., in
MIMO In radio, multiple-input and multiple-output, or MIMO (), is a method for multiplying the capacity of a radio link using multiple transmission and receiving antennas to exploit multipath propagation. MIMO has become an essential element of wi ...
systems) or cyclo-stationary signals.


See also

*
Split Levinson recursion Split(s) or The Split may refer to: Places * Split, Croatia, the largest coastal city in Croatia * Split Island, Canada, an island in the Hudson Bay * Split Island, Falkland Islands * Split Island, Fiji, better known as Hạfliua Arts, entertai ...
* Linear prediction *
Autoregressive model In statistics, econometrics and signal processing, an autoregressive (AR) model is a representation of a type of random process; as such, it is used to describe certain time-varying processes in nature, economics, etc. The autoregressive model spe ...


Notes


References

Defining sources * Levinson, N. (1947). "The Wiener RMS error criterion in filter design and prediction." ''J. Math. Phys.'', v. 25, pp. 261–278. * Durbin, J. (1960). "The fitting of time series models." ''Rev. Inst. Int. Stat.'', v. 28, pp. 233–243. * Trench, W. F. (1964). "An algorithm for the inversion of finite Toeplitz matrices." ''J. Soc. Indust. Appl. Math.'', v. 12, pp. 515–522. * Musicus, B. R. (1988). "Levinson and Fast Choleski Algorithms for Toeplitz and Almost Toeplitz Matrices." ''RLE TR'' No. 538, MIT

* Delsarte, P. and Genin, Y. V. (1986). "The split Levinson algorithm." ''IEEE Transactions on Acoustics, Speech, and Signal Processing'', v. ASSP-34(3), pp. 470–478. Further work * * Richard P. Brent, Brent R.P. (1999), "Stability of fast algorithms for structured linear systems", ''Fast Reliable Algorithms for Matrices with Structure'' (editors—T. Kailath, A.H. Sayed), ch.4 (
SIAM Thailand ( ), historically known as Siam () and officially the Kingdom of Thailand, is a country in Southeast Asia, located at the centre of the Indochinese Peninsula, spanning , with a population of almost 70 million. The country is bo ...
). * Bunch, J. R. (1985). "Stability of methods for solving Toeplitz systems of equations." ''SIAM J. Sci. Stat. Comput.'', v. 6, pp. 349–364

* Summaries * Bäckström, T. (2004). "2.2. Levinson–Durbin Recursion." ''Linear Predictive Modelling of Speech – Constraints and Line Spectrum Pair Decomposition.'' Doctoral thesis. Report no. 71 / Helsinki University of Technology, Laboratory of Acoustics and Audio Signal Processing. Espoo, Finland

* Claerbout, Jon F. (1976). "Chapter 7 – Waveform Applications of Least-Squares." ''Fundamentals of Geophysical Data Processing.'' Palo Alto: Blackwell Scientific Publications

*{{Citation , last1=Press, first1=WH, last2=Teukolsky, first2=SA, last3=Vetterling, first3=WT, last4=Flannery, first4=BP, year=2007, title=Numerical Recipes: The Art of Scientific Computing, edition=3rd, publisher=Cambridge University Press, publication-place=New York, isbn=978-0-521-88068-8, chapter=Section 2.8.2. Toeplitz Matrices, chapter-url=http://apps.nrbook.com/empanel/index.html?pg=96 * Golub, G.H., and Loan, C.F. Van (1996). "Section 4.7 : Toeplitz and related Systems" ''Matrix Computations'', Johns Hopkins University Press Matrices Numerical analysis