Levinson-Durbin algorithm

• We pursue a block decomposition of $R_{A,p}$. Subsitituting into the inversion formula (16)
  

  $$A = R_{A,p-1}$$

  $$B = C^{*} = \rho_{p-1}$$

  $$D = R(0)$$

  
  obtains:
  

  $$R^{-1}_{A,p} = \left[\begin{array}{cc} R_{A,p-1} & \rho^{*\char93 }_{p-1} \\ \rho^{\char93 }_{p-1} & R(0) \end{array}\right]^{-1}$$

  $$= \begin{normalsize} \left[\begin{array}{cc} R^{-1}_{A,p-1} + R^{-1... ... }_{p-1}R^{-1}_{A,p-1} & \lambda^{-1}_{p-1} \end{array}\right] \end{normalsize}$$

  
  

• The Schur complement, denoted as $\lambda_{p-1}$, is nothing but the PSD-minimal squared error for the model of order $p-1$:
  

  $$\lambda_{p-1} = R(0) - \rho^{\char93 }_{p-1}R^{-1}_{A,p-1}\rho^{*\char93 }_{p-1}$$

  $$= R(0) + A^{\char93 }_{p-1}\rho^{*\char93 }_{p-1}$$

  $$= R(0) + A_{1,p-1}R(1) + \ldots + A_{p-1,p-1}R(p-1)$$

  $$= \sum y_{n}\left[y_{n} + A_{1,p-1}y_{n-1} + \ldots + A_{p-1,p-1}y_{n-(p-1)}\right]^{*}$$

  $$= \sum y_{n}e^{*}_{n,p-1}$$

  $$= \sum e_{n,p-1}e^{*}_{n,p-1}$$ (19)

  
  The last step follows by writing $y_{n} = e_{n,p-1}-A_{1,p-1}y_{n-1}- \ldots A_{p-1,p-1}y_{n-(p-1)}$ and using orthogonality of the error $e_{n,p-1}$ and past outputs $y_{n-1} \ldots y_{n-(p-1)}$ under optimal choice of model parameters, i.e.
  

  $$\sum y_{n-k}e^{*}_{n,p-1} = {\bf0}, \quad k = 1 \ldots p-1$$ (20)

  
  See Appendix for proof. Note that $\lambda_{p}/N$ gives a (biased) estimate of the prediction error covariance of order $p$.

• From the formula for $R^{-1}_{A,p}$, we find a recursion for $A_{p}$:
  

  $$A_{p} = \rho_{p} R^{-1}_{A,p}$$

  $$= \begin{normalsize} \left[\begin{array}{c\vert c} \rho_{p-1} & R(p... ... }_{p-1}R^{-1}_{A,p-1} & \lambda^{-1}_{p-1} \end{array}\right] \end{normalsize}$$

  $$= \begin{normalsize} \left[\begin{array}{c\vert c} \rho_{p-1}R^{-1}... ...{*\char93 }_{p-1}\right) \lambda^{-1}_{p-1} \end{array}\right] \end{normalsize}$$

  $$= \left[\begin{array}{c\vert c} A_{p-1} - k_{p}A^{\char93 }_{p-1} & k_{p}\end{array}\right]$$ (21)

  
  In the final step, we define $k_{p}$ = $\Delta_{p-1} \lambda^{-1}_{p-1}$, where $\Delta_{p-1} = R(p) - \rho^{*\char93 }_{p-1}R^{-1}_{p-1} \rho^{\char93 }_{p-1}$, like $\lambda_{p-1}$, is in the form of a Schur complement.

• To interpret $\Delta_{p-1}$ and $k_{p}$, we parallel the development for $\lambda_{p-1}$ (19):
  

  $$\Delta_{p-1} = R(p) - \rho_{p-1}R^{-1}_{A,p-1}\rho^{*\char93 }_{p-1}$$

  $$= R(p) + A_{1,p-1}R(p-1) + \ldots + A_{p-1,p-1}R(1)$$

  $$= \sum y_{n}\left[y_{n-p} + A_{1,p-1}y_{n-(p-1)} + \ldots + A_{p-1,p-1}y_{n-1)}\right]^{*}$$

  
  The quantity $y_{n-p} + A_{1,p-1}y_{n-(p-1)} + \ldots + A_{p-1,p-1}y_{n-1}$ is interpreted as a backwards prediction error, and we denote as $f_{n,p-1}$. Proceeding:
  

  $$\Delta_{p-1} = \sum y_{n}f^{*}_{n,p-1}$$

  $$= \sum e_{n,p-1}f^{*}_{n,p-1}$$ (22)

  
  Again by orthogonality arguments, (see Appendix), we replace $y_{n}$ by $e_{n,p-1}$. Note that $Delta_{p}/N$ gives a (biased) estimate of the cross-covariance between forward and backward prediction errors of order $p$.

• The recursion for $A_{p}$ depends on the unknown quantities $\lambda_{p-1}$ and $\Delta_{p-1}$. Because $\lambda_{p-1} = R(0) - A_{p-1}\rho^{*}_{p-1}$ and $\Delta_{p-1} = R(p) - A_{p-1}\rho^{*\char93 }_{p-1}$, we see everything is available for the computation. However, forming $A_{p-1}\rho^{*}_{p-1}$ takes $p-1$ $m$-by-$m$ matrix multiplications. To eliminate all but one of these multiplications, we develop a recursion for $\lambda_{p}$:
  

  $$\lambda_{p} = R(0) - A_{p}\rho^{*}_{p}$$

  $$= R(0) - \left[\begin{array}{c\vert c} A_{p-1} - k_{p}A^{\char93 }_... ...}\end{array}\right] \left[\begin{array}{c} \rho_{p-1} \\ R(p)\end{array}\right]$$

  $$= \left(R(0) - A_{p-1}\rho_{p-1}\right) - k_{p}\left(R(0) - A_{p-1}\rho_{p-1}\right)$$

  $$= \lambda_{p-1} - k_{p}\Delta^{*}_{p-1}$$ (23)

  
  This gives the standard form of the Levinson-Durbin algorithm.

• For the initialization, consider $p=0$. Since no prediction is done, the forward error is nothing but the observation: $e_{n,0} = y_{n}$ and the backward error is nothing but the past observation: $f_{n,0} = y_{n-1}$, given the way it was defined (23). It follows that $\lambda_{0} = R(0)$ and $\Delta_{0}$ = $R(1)$; consequently, $k_{1} = R(1)\left[R(0)\right]^{-1}$, which equals $A_{1}$ thanks to the original least squares solution (12). The Levinson-Durbin algorithm may be summarized:
  

  $$\lambda_{p} = \lambda_{p-1} - k_{p}\Delta^{*}_{p-1}$$

  $$\Delta_{p} = R(p) - A_{p-1}\rho^{*\char93 }_{p-1}$$

  $$k_{p+1} = \Delta_{p} \lambda^{-1}_{p}$$

  $$A_{p+1} = \left[\begin{array}{c\vert c} A_{p} - k_{p+1}A^{\char93 }_{p} & k_{p+1}\end{array}\right]$$ Initialization:

  $$\lambda_{0} = R(0)$$

  $$\Delta_{0} = R(1)$$

  $$k_{1} = \Delta_{0} \lambda^{-1}_{0}$$

  $$A_{1} = k_{1}$$