Levinson recursion
|
|
Levinson recursion is a mathematical procedure which recursively calculates the solution to a Toeplitz matrix.
It was proposed by Norman Levinson in 1947, improved by Durbin in 1960, and given a matrix formulation requiring <math>4n^2<math> multiplications by W. F. Trench in 1964. S. Zohar improved Trench's algorithm in 1969, requiring only <math>3n^2<math> multiplications.
In most applications where Toeplitz matrices appear, we are dealing with a problem such as <math>\mathbf{Ta=b}<math> where <math>\mathbf T<math> is an <math>n\times n<math> Toeplitz matrix and <math>\mathbf a<math> and <math>\mathbf b<math> are vectors. The problem is to solve <math>\mathbf a<math> when <math>\mathbf T<math> and <math>\mathbf b<math> are known. For the solution, straightforward application of Gaussian elimination is rather inefficient, with complexity, <math>O(n^3)<math>, since it does not employ the strong structures present in the Toeplitz system.
A first improvement to Gaussian elimination is the Levinson recursion which can be applied to symmetric Toeplitz systems. To illustrate the basics of the Levinson algorithm, first define the <math>p\times p<math> principal sub-matrix <math>T_p<math> as the upper left block of <math>\mathbf T<math>. Further, assume that we have the order <math>p<math> solution <math>{\mathbf a}_p<math> to equation
- <math>
\begin{bmatrix}
t_0 & t_1 & t_2 & \dots & t_p \\
t_{1} & t_0 & t_1 & \ddots & t_{p-1} \\
t_{2} & t_{1} & t_0 & \ddots & t_{p-2} \\
\vdots & \ddots & \ddots & \ddots & \vdots \\
t_{p} & t_{p-1} & t_{p-2} & \dots & t_{0} \\
\end{bmatrix}
\begin{bmatrix}
1 \\ a^{(p)}_1 \\ \vdots \\ a^{(p)}_p
\end{bmatrix}
=
\begin{bmatrix}
\epsilon_p \\ 0 \\ \vdots \\ 0
\end{bmatrix}.
<math> Extension of <math>{\mathbf a}_p<math> with a zero yields
- <math>
\begin{bmatrix}
t_0 & t_1 & t_2 & \dots & t_{p+1} \\
t_{1} & t_0 & t_1 & \ddots & t_{p } \\
t_{2} & t_{1} & t_0 & \ddots & t_{p-1} \\
\vdots & \ddots & \ddots & \ddots & \vdots \\
t_{p+1}& t_{p } & t_{p-1} & \dots & t_{0} \\
\end{bmatrix}
\begin{bmatrix}
1 \\ a^{(p)}_1 \\ \vdots \\ a^{(p)}_p \\ 0
\end{bmatrix}
=
\begin{bmatrix}
\epsilon_p \\ 0 \\ \vdots \\ 0 \\ \eta_p
\end{bmatrix}
<math> where <math>\eta_p=\sum_{i=0}^pa_i^{(p)}t_{p-i+1}<math> and <math>a_0^{(p)}=1<math>. The salient step comes through the realisation that since <math>{\mathbf T}_p{\mathbf a}_p={\mathbf u}_p<math> where
<math>{\mathbf u}_p=[\epsilon_p\,0\,\dots\,0]^T<math> and <math>{\mathbf T}_p<math> is
symmetric, we have <math>{\mathbf T}_p{\mathbf a}_p^\#={\mathbf u}_p^\#<math>, where superscript # denotes reversal of rows. By defining a reflection coefficient <math>\Gamma_p<math>, we obtain
- <math>
{\mathbf T}_{p+1}
\left(
\begin{bmatrix}
1 \\ a^{(p)}_1 \\ \vdots \\ a^{(p)}_p \\ 0
\end{bmatrix}
+
\Gamma_p
\begin{bmatrix}
0 \\ a^{(p)}_p \\ a^{(p)}_{p-1} \\ \vdots \\ 1
\end{bmatrix}
\right)
=
\left(
\begin{bmatrix}
\epsilon_p \\ 0 \\ \vdots \\ 0 \\ \eta_p
\end{bmatrix}
+\Gamma_p
\begin{bmatrix}
\eta_p \\ 0 \\ \vdots \\ 0 \\ \epsilon_p
\end{bmatrix}
\right).<math>
Obviously, choosing <math>\Gamma_p<math> so that <math>\eta_p+\Gamma_p\epsilon_p=0<math> yields the order <math>p+1<math> solution to <math>{\mathbf{Ta=b}}<math> as
- <math>
{\mathbf a}_{p+1}=\begin{bmatrix}{\mathbf a}_p\\0\end{bmatrix} +\Gamma_p\begin{bmatrix}0\\{\mathbf a}_p^\#\end{bmatrix}. <math> Consequently, with a suitable choice of initial values (<math>{\mathbf a}_0=1<math>), this procedure can be used to recursively solve equations of type <math>\mathbf{Ta}=[\sigma^2\,0\,0\dots0]^T<math>. Furthermore, using the intermediate values <math>{\mathbf a}_p<math>, it is straightforward to solve arbitrary problems of type <math>\mathbf{Ta=b}<math>. The former algorithm, is often called the Levinson-Durbin recursion and the latter, the solution of arbitrary Toeplitz equations, the Levinson recursion.
While the Levinson-Durbin recursion has the complexity of <math>O(n^2)<math>, it is possible to further improve the algorithm to reduce complexity by half. The algorithm, called the split Levinson-Durbin algorithm, uses a three term recursion instead of the two term recursion in conventional Levinson recursion. Then either the symmetric or antisymmetric part of two consecutive order solutions, <math>{\mathbf a}_{p-1}<math> and <math>{\mathbf a}_{p}<math> are used to obtain the next order solution <math>{\mathbf a}_{p+1}<math>.
Several other possibilities exist for the solution of Toeplitz systems, such as the Schur or Cholesky decompositions, but the split Levinson-Durbin is the most efficient. The others are sometimes preferred when decimal truncations cause numerical instability.
References
- Bunch, J. R. (1985). "Stability of methods for solving Toeplitz systems of equations." SIAM J. Sci. Stat. Comput., v. 6, pp. 349-364.
- Trench, W. F. (1964). "An algorithm for the inversion of finite Toeplitz matrices." J. Soc. Indust. Appl. Math., v. 12, pp. 515-522.