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

4n^2

multiplications by W. F. Trench in 1964. S. Zohar improved Trench's algorithm in 1969, requiring only

3n^2

multiplications. In most applications where Toeplitz matrices appear, we are dealing with a problem such as

\mathbf{Ta=b}

where

\mathbf T

is an

n\times n

Toeplitz matrix and

\mathbf a

and

\mathbf b

are vectors. The problem is to solve

\mathbf a

when

\mathbf T

and

\mathbf b

are known. For the solution, straightforward application of Gaussian elimination is rather inefficient, with complexity,

O(n^3)

, 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

p\times p

principal sub-matrix

T_p

as the upper left block of

\mathbf T

. Further, assume that we have the order

p

solution

{\mathbf a}_p

to equation

   \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}.

Extension of

{\mathbf a}_p

with a zero yields

   \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}

where

\eta_p=\sum_{i=0}^pa_i^{(p)}t_{p-i+1}

and

a_0^{(p)}=1

. The salient step comes through the realisation that since

{\mathbf T}_p{\mathbf a}_p={\mathbf u}_p

where

   ${\mathbf u}_p= \epsilon_p\,0\,\dots\,0^T$  and  ${\mathbf T}_p$  is

symmetric, we have

{\mathbf T}_p{\mathbf a}_p^\#={\mathbf u}_p^\#

, where superscript # denotes reversal of rows. By defining a reflection coefficient

\Gamma_p

, we obtain

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

Obviously, choosing

\Gamma_p

so that

\eta_p+\Gamma_p\epsilon_p=0

yields the order

p+1

solution to

{\mathbf{Ta=b}}

{\mathbf a}_{p+1}=\begin{bmatrix}{\mathbf a}_p\\0\end{bmatrix} +\Gamma_p\begin{bmatrix}0\\{\mathbf a}_p^\#\end{bmatrix}. Consequently, with a suitable choice of initial values (

{\mathbf a}_0=1

), this procedure can be used to recursively solve equations of type

\mathbf{Ta}= \sigma^2\,0\,0\dots0^T

. Furthermore, using the intermediate values

{\mathbf a}_p

, it is straightforward to solve arbitrary problems of type

\mathbf{Ta=b}

. 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

O(n^2)

, 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,

{\mathbf a}_{p-1}

and

{\mathbf a}_{p}

are used to obtain the next order solution

{\mathbf a}_{p+1}

. 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.

Levinson Recursion

References

See also