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 [itex]4n^2[itex] multiplications by W. F. Trench in 1964. S. Zohar improved Trench's algorithm in 1969, requiring only [itex]3n^2[itex] multiplications.

In most applications where Toeplitz matrices appear, we are dealing with a problem such as [itex]\mathbf{Ta=b}[itex] where [itex]\mathbf T[itex] is an [itex]n\times n[itex] Toeplitz matrix and [itex]\mathbf a[itex] and [itex]\mathbf b[itex] are vectors. The problem is to solve [itex]\mathbf a[itex] when [itex]\mathbf T[itex] and [itex]\mathbf b[itex] are known. For the solution, straightforward application of Gaussian elimination is rather inefficient, with complexity, [itex]O(n^3)[itex], 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 [itex]p\times p[itex] principal sub-matrix [itex]T_p[itex] as the upper left block of [itex]\mathbf T[itex]. Further, assume that we have the order [itex]p[itex] solution [itex]{\mathbf a}_p[itex] to equation

[itex]
\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}.

[itex] Extension of [itex]{\mathbf a}_p[itex] with a zero yields

[itex]
\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}

[itex] where [itex]\eta_p=\sum_{i=0}^pa_i^{(p)}t_{p-i+1}[itex] and [itex]a_0^{(p)}=1[itex]. The salient step comes through the realisation that since [itex]{\mathbf T}_p{\mathbf a}_p={\mathbf u}_p[itex] where

[itex]{\mathbf u}_p=[\epsilon_p\,0\,\dots\,0]^T[itex] and [itex]{\mathbf T}_p[itex] is

symmetric, we have [itex]{\mathbf T}_p{\mathbf a}_p^\#={\mathbf u}_p^\#[itex], where superscript # denotes reversal of rows. By defining a reflection coefficient [itex]\Gamma_p[itex], we obtain

[itex]
{\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).[itex]

Obviously, choosing [itex]\Gamma_p[itex] so that [itex]\eta_p+\Gamma_p\epsilon_p=0[itex] yields the order [itex]p+1[itex] solution to [itex]{\mathbf{Ta=b}}[itex] as

[itex]

{\mathbf a}_{p+1}=\begin{bmatrix}{\mathbf a}_p\\0\end{bmatrix} +\Gamma_p\begin{bmatrix}0\\{\mathbf a}_p^\#\end{bmatrix}. [itex] Consequently, with a suitable choice of initial values ([itex]{\mathbf a}_0=1[itex]), this procedure can be used to recursively solve equations of type [itex]\mathbf{Ta}=[\sigma^2\,0\,0\dots0]^T[itex]. Furthermore, using the intermediate values [itex]{\mathbf a}_p[itex], it is straightforward to solve arbitrary problems of type [itex]\mathbf{Ta=b}[itex]. 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 [itex]O(n^2)[itex], 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, [itex]{\mathbf a}_{p-1}[itex] and [itex]{\mathbf a}_{p}[itex] are used to obtain the next order solution [itex]{\mathbf a}_{p+1}[itex].

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.

• Art and Cultures
• Countries of the World (http://www.academickids.com/encyclopedia/index.php/Countries)
• Space and Astronomy