Editing Levinson recursion (section)

== Derivation ==

=== Background ===
Matrix equations follow the form

: <math>\mathbf M \, \vec x = \vec y.</math>

The Levinson–Durbin algorithm may be used for any such equation, as long as '''M''' is a known [[Toeplitz matrix]] with a nonzero main diagonal. Here <math>\vec y</math> is a known [[vector space|vector]], and <math>\vec x</math> is an unknown vector of numbers ''x''<sub>''i''</sub> yet to be determined.

For the sake of this article, ''ê''<sub>''i''</sub> 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'''<sup>''n''</sup> is an ''n''×''n'' matrix that copies the upper left ''n''×''n'' block from '''M''' – that is, ''T''<sup>''n''</sup><sub>''ij''</sub> = ''M''<sub>''ij''</sub>.

'''T'''<sup>''n''</sup> is also a Toeplitz matrix, meaning that it can be written as

: <math>\mathbf T^n = \begin{bmatrix}
    t_0    & t_{-1}  & t_{-2}  & \dots  & t_{-n+1}   \\
    t_1    & t_0     & t_{-1}  & \dots  & t_{-n+2} \\
    t_2    & t_1     & t_0     & \dots  & t_{-n+3} \\
    \vdots & \vdots  & \vdots  & \ddots & \vdots   \\
    t_{n-1}& t_{n-2} & t_{n-3} & \dots  & t_0
  \end{bmatrix}. </math>

=== 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''<sup>th</sup> "forward vector", denoted <math>\vec f^n</math>, as the vector of length ''n'' which satisfies:

:<math>\mathbf T^n \vec f^n = \hat e_1.</math>

The ''n''<sup>th</sup> "backward vector" <math>\vec b^n</math> is defined similarly; it is the vector of length ''n'' which satisfies:

:<math>\mathbf T^n \vec b^n = \hat e_n.</math>

An important simplification can occur when '''''M''''' is a [[symmetric matrix]]; then the two vectors are related by ''b''<sup>''n''</sup><sub>''i''</sub> = ''f''<sup>''n''</sup><sub>''n''+1−''i''</sub>—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''<sup>th</sup> forward and backward vector may be found from the vectors of length ''n''&nbsp;−&nbsp;1 as follows. First, the forward vector may be extended with a zero to obtain:

:<math>\mathbf T^n \begin{bmatrix} \vec f^{n-1} \\ 0 \\ \end{bmatrix} =
  \begin{bmatrix}
    \        & \               & \     & t_{-n+1}   \\
    \        & \mathbf T^{n-1} & \     & t_{-n+2} \\
    \        & \               & \     & \vdots   \\
    t_{n-1}  & t_{n-2}         & \dots & t_0      \\
  \end{bmatrix}
  \begin{bmatrix}  \            \\
                   \vec f^{n-1} \\
                   \            \\
                   0            \\
                   \            \\
  \end{bmatrix} = 
  \begin{bmatrix}  1            \\
                   0            \\
                   \vdots       \\
                   0            \\
                   \varepsilon_f^n
  \end{bmatrix}. </math>

In going from '''''T''<sup>''n''−1</sup>''' to '''''T<sup>n</sup>''''', 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 ''ε''<sub>''f''</sub> which occurs in the last place. The above equation gives it the value of:

: <math> \varepsilon_f^n \ = \  \sum_{i=1}^{n-1} \ M_{ni} \  f_{i}^{n-1} \ = \ \sum_{i=1}^{n-1} \  t_{n-i} \ f_{i}^{n-1}. </math>

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,

:<math> \mathbf T^n \begin{bmatrix} 0 \\ \vec b^{n-1} \\ \end{bmatrix} =
\begin{bmatrix}
    t_0     & \dots & t_{-n+2}         & t_{-n+1} \\
    \vdots  & \     & \               & \       \\
    t_{n-2} & \     & \mathbf T^{n-1} & \       \\
    t_{n-1} & \     & \               &
  \end{bmatrix}
  \begin{bmatrix}  \            \\
                   0            \\
                   \            \\
                   \vec b^{n-1} \\
                   \            \\
  \end{bmatrix} = 
  \begin{bmatrix}  \varepsilon_b^n  \\
                   0             \\
                   \vdots        \\
                   0             \\
                   1
  \end{bmatrix}. </math>

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 ''ε''<sub>''b''</sub> with value:

:<math> \varepsilon_b^n \ = \ \sum_{i=2}^n \  M_{1i} \ b_{i-1}^{n-1} \ = \ \sum_{i=1}^{n-1} \  t_{-i} \ b_i^{n-1}. \ </math>

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 <math>(\alpha,\beta)</math>:

:<math>\mathbf T \left( \alpha
  \begin{bmatrix}
                   \vec f \\
                   \            \\
                   0            \\
  \end{bmatrix} + \beta
  \begin{bmatrix}
                   0            \\
                   \            \\
                   \vec b
  \end{bmatrix} \right ) = \alpha
  \begin{bmatrix}  1        \\
                   0        \\
                   \vdots   \\
                   0        \\
                   \varepsilon_f \\
  \end{bmatrix} + \beta
  \begin{bmatrix}  \varepsilon_b  \\
                   0             \\
                   \vdots        \\
                   0             \\
                   1
  \end{bmatrix}.</math>

If ''α'' and ''β'' are chosen so that the right hand side yields ''ê''<sub>1</sub> or ''ê''<sub>''n''</sub>, then the quantity in the parentheses will fulfill the definition of the ''n''<sup>th</sup> 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, <math>\alpha^n_{f}</math>, <math>\beta^n_{f}</math> are such that :
:<math>
\vec f^n = \alpha^n_{f} \begin{bmatrix} \vec f^{n-1}\\
0
\end{bmatrix}
+\beta^n_{f}\begin{bmatrix}0\\
\vec b^{n-1}
\end{bmatrix}
</math>
and respectively  <math>\alpha^n_{b}</math>, <math>\beta^n_{b}</math> are such that :
:<math>\vec b^n = \alpha^n_{b}
\begin{bmatrix}
\vec f^{n-1}\\
0
\end{bmatrix}
+\beta^n_{b}\begin{bmatrix}
0\\
\vec b^{n-1}
\end{bmatrix}.
</math>
By multiplying both previous equations by <math>{\mathbf T}^n</math> one gets the following equation:
: <math>
\begin{bmatrix} 1 & \varepsilon^n_b \\
0 & 0 \\
\vdots & \vdots \\
0 & 0 \\
\varepsilon^n_f & 1
\end{bmatrix} \begin{bmatrix} \alpha^n_f & \alpha^n_b \\ \beta^n_f & \beta^n_b \end{bmatrix}
= \begin{bmatrix}
1 & 0 \\
0 & 0 \\
\vdots & \vdots \\
0 & 0 \\
0 & 1
\end{bmatrix}.</math>

Now, all the zeroes in the middle of the two vectors above being disregarded and collapsed, only the following equation is left:

: <math> \begin{bmatrix} 1 & \varepsilon^n_b \\ \varepsilon^n_f & 1 \end{bmatrix} \begin{bmatrix} \alpha^n_f & \alpha^n_b \\ \beta^n_f & \beta^n_b \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}.</math>

With these solved for (by using the Cramer 2×2 matrix inverse formula), the new forward and backward vectors are:

: <math>\vec f^n = {1 \over { 1 - \varepsilon_b^n \varepsilon_f^n }}          \begin{bmatrix} \vec f^{n-1} \\ 0 \end{bmatrix}
                 - { \varepsilon_f^n \over { 1 - \varepsilon_b^n \varepsilon_f^n }}\begin{bmatrix} 0 \\ \vec b^{n-1} \end{bmatrix}</math>

: <math>\vec b^n =  { 1 \over { 1 - \varepsilon_b^n \varepsilon_f^n }}         \begin{bmatrix} 0 \\ \vec b^{n-1} \end{bmatrix}
                 - { \varepsilon_b^n \over { 1 - \varepsilon_b^n \varepsilon_f^n }}\begin{bmatrix} \vec f^{n-1} \\ 0 \end{bmatrix}.</math>

Performing these vector summations, then, gives the ''n''<sup>th</sup> 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:

: <math>\vec f^1 = \vec b^1 = \left[ {1 \over M_{11}} \right] = \left[ {1 \over t_0} \right].</math>

=== Using the backward vectors ===
The above steps give the ''N'' backward vectors for '''''M'''''. From there, a more arbitrary equation is:

: <math> \vec y = \mathbf M \  \vec x. </math>

The solution can be built in the same recursive way that the backwards vectors were built. Accordingly, <math>\vec x</math> must be generalized to a sequence of intermediates <math>\vec x^n</math>, such that <math>\vec x^N = \vec x</math>.

The solution is then built recursively by noticing that if

: <math> \mathbf T^{n-1}
  \begin{bmatrix}  x_1^{n-1}     \\
                   x_2^{n-1}     \\
                   \vdots   \\
                   x_{n-1}^{n-1} \\
  \end{bmatrix} =   
  \begin{bmatrix}  y_1     \\
                   y_2     \\
                   \vdots   \\
                   y_{n-1}
  \end{bmatrix}.</math>

Then, extending with a zero again, and defining an error constant where necessary:

: <math> \mathbf T^{n}
  \begin{bmatrix}  x_1^{n-1}     \\
                   x_2^{n-1}     \\
                   \vdots   \\
                   x_{n-1}^{n-1} \\
                   0
  \end{bmatrix} =   
  \begin{bmatrix}  y_1     \\
                   y_2     \\
                   \vdots   \\
                   y_{n-1} \\
                   \varepsilon_x^{n-1}
  \end{bmatrix}.</math>

We can then use the ''n''<sup>th</sup> backward vector to eliminate the error term and replace it with the desired formula as follows:

: <math> \mathbf T^{n} \left (
  \begin{bmatrix}  x_1^{n-1}     \\
                   x_2^{n-1}     \\
                   \vdots   \\
                   x_{n-1}^{n-1} \\
                   0 \\
  \end{bmatrix} + (y_n - \varepsilon_x^{n-1}) \  \vec b^n \right ) =   
  \begin{bmatrix}  y_1     \\
                   y_2     \\
                   \vdots   \\
                   y_{n-1} \\
                   y_n
  \end{bmatrix}.</math>

Extending this method until ''n'' = ''N'' yields the solution <math>\vec x</math>.

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.