# Algoritmo de Levinson

El algoritmo de Levinson o de Levinson-Durbin es un algoritmo del álgebra lineal para calcular en forma recursiva la solución de una ecuación que involucra una matriz de Toeplitz. El costo computacional es de Θ(n$^2$), una mejora considerable frente a la eliminación de Gauss-Jordan, cuyo costo es de Θ(n$^3$).
Hay algoritmos nuevos, llamados asintóticamente rápidos o a menudo algoritmos de Toeplitz superrápidos, que pueden resolver con un costo de Θ(n log$^p$n) para varios p (por ejemplo, para p = 2, p = 3). La recursión de Levinson sigue siendo popular por distintas razones; por un lado, es relativamente simple de comprender en comparación; por otro lado, puede ser más rápida que un algoritmo superrápido para n pequeño (normalmente para n < 256 Archivado el 5 de septiembre de 2006 en Wayback Machine.).
El algoritmo de Levinson-Durbin fue propuesto por primera vez por Norman Levinson en 1947, mejorado por J. Durbin en 1960 y más tarde mejorado a 4n$^2$ y luego a 3n$^2$ multiplicaciones por W. F. Trench y S. Zohar, respectivamente.
Otros métodos para procesar datos incluyen la descomposición de Schur y la descomposición de Cholesky. En comparación a estos, la recursión de Levinson (particularmente la recursión de Split-Levinson) tiende a ser más rápida computacionalmente, aunque más sensible a imperfecciones computacionales como errores de redondeo.

## Derivación
### Contexto

Las ecuaciones matriciales tienen la siguiente forma:
    ${\vec {y}}=\mathbf {M} \ {\vec {x}}$
El algoritmo de Levinson-Durbin puede usarse para tales ecuaciones, siempre y cuando M sea una matriz de Toeplitz conocida, con una diagonal principal no nula. En la ecuación, ${\vec {y}}$ es un vector conocido (en el sentido del álgebra lineal de lista de números) y ${\vec {x}}$ es un vector de valores x$_i$ desconocidos y a ser determinados.
Llamaremos ê$_i$ al vector compuesto de ceros, salvo en su i-ésimo lugar, que contiene el valor de la unidad. Su dimensión dependerá del contexto y estará implícita en éste. Con N nos referiremos al ancho de la matriz anterior, es decir, M es una matriz de N×N. Finalmente, los superindices se referirán a índices inductivos, mientras que subíndices denotarán simplemente índices. Por ejemplo, la matriz T$^n$ es una matriz de n×n que copia el bloque superior izquierdo de n×n de M, es decir, T$^n$$_{ij}$ = M$_{ij}$.
T$^n$ también es una matriz de Toeplitz, en el sentido de que puede escribirse como:
    $\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}}$

### Pasos introductorios

El algoritmo se desarrolla en dos pasos. En el primero, se establecen dos conjuntos de vectores, llamados vectores de avance y de retroceso. Los vectores de avance sirven para obtener el conjunto de vectores de retroceso y pueden luego ser descartados. Los vectores de retroceso son necesarios para el segundo paso, donde se los usa para construir la solución.
La recursión de Levinson-Durbin define el n-ésimo "vector de avance", denotado ${\vec {f}}^{n}$, como el vector de longitud n que satisface:
    $\mathbf {T} ^{n}{\vec {f}}^{n}={\hat {e}}_{1}$
El n-ésimo "vector de retroceso" ${\vec {b}}^{n}$ se define de manera similar; es el vector de longitud n que satisface:
    $\mathbf {T} ^{n}{\vec {b}}^{n}={\hat {e}}_{n}$
Una simplificación importante ocurre cuando M es una matriz simétrica: los dos vectores se relacionan mediante b$^n$$_i$ = f$^n$$_{n+1-i}$; es decir, uno se obtiene invirtiendo el order de los elementos del otro. Esto puede ahorrar cálculos extra en ese caso en particular.

### Obtención del vector de retroceso

Incluso si la matriz no es simétrica, los n-ésimos vectores de avance y retroceso pueden encontrarse a partir de los vectores de longitud n-1 de la siguiente manera. Primero, el vector de avance puede extenderse con un cero para obtener:
    $\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\\\epsilon _{f}^{n}\\\end{bmatrix}}$
Al ir de T$^{n-1}$ a T$^n$, la columna extra agregada a la matriz no altera el resultado cuando se usa un cero para extender el vector de avance. Sin embargo, la fila extra agregada a la matriz sí lo altera, creando un error indeseado ε$_f$ en el último lugar. De la ecuación anterior:
    $\epsilon _{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}$
Volveremos a este error en un momento. Antes, extenderemos al vector de retroceso de manera similar (aunque invertida). Para el vector de retroceso:
    $\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}\epsilon _{b}^{n}\\0\\\vdots \\0\\1\\\end{bmatrix}}$
Como antes, la columna extra agregada a la matriz no altera el resultado, pero sí lo hace la fila extra. Tenemos aquí, entonces, otro error indeseado ε$_b$, cuyo valor es:
    $\epsilon _{b}^{n}\ =\ \sum _{i=2}^{n}\ M_{1i}\ b_{i-1}^{n-1}\ =\ \sum _{i=1}^{n-1}\ t_{-i}\ b_{i}^{n-1}\$
Podemos usar estos dos términos de error para eliminarlos en forma simultánea. Debido a la propiedad de linealidad de las matrices:
    $\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\\\epsilon _{f}\\\end{bmatrix}}+\beta {\begin{bmatrix}\epsilon _{b}\\0\\\vdots \\0\\1\\\end{bmatrix}}$
Si elegimos α y β de manera tal que el lado derecho dé como resultado ê$_1$ o ê$_n$, entonces el valor entre paréntesis va a satisfacer la definición del vector de avance o retroceso n-ésimo, respectivamente. Con esos valores de α y β, la suma de vectores entre paréntesis es simple y da el resultado esperado.
Para encontrar estos coeficientes, es decir, α$_f$, β$_f$, α$_b$, y β$_b$, nótese que todos los ceros en el medio de los dos vectores de arriba pueden ser ignorados y colapsados, dejando simplemente la ecuación:
    ${\begin{bmatrix}1&\epsilon _{b}\\\epsilon _{f}&1\end{bmatrix}}{\begin{bmatrix}\alpha _{f}&\alpha _{b}\\\beta _{f}&\beta _{b}\end{bmatrix}}={\begin{bmatrix}1&0\\0&1\end{bmatrix}}$
Resolviendo este sistema (por ejemplo, con la fórmula de inversa para matrices de 2x2), los nuevos vectores de avance y retroceso son:
    ${\vec {f}}^{n}={1 \over {1-\epsilon _{b}^{n}\epsilon _{f}^{n}}}{\begin{bmatrix}{\vec {f}}^{n-1}\\0\end{bmatrix}}-{\epsilon _{f}^{n} \over {1-\epsilon _{b}^{n}\epsilon _{f}^{n}}}{\begin{bmatrix}0\\{\vec {b}}^{n-1}\end{bmatrix}}$
    ${\vec {b}}^{n}={1 \over {1-\epsilon _{b}^{n}\epsilon _{f}^{n}}}{\begin{bmatrix}0\\{\vec {b}}^{n-1}\end{bmatrix}}-{\epsilon _{b}^{n} \over {1-\epsilon _{b}^{n}\epsilon _{f}^{n}}}{\begin{bmatrix}{\vec {f}}^{n-1}\\0\end{bmatrix}}$
Realizando la suma de vectores, entonces, se obtienen los vectores de avance y retroceso n-ésimos, dados los anteriores. Lo único que resta es encontrar los vectores iniciales, para luego obtener los otros mediante simples sumas y multiplicaciones. Los vectores de avance y retroceso iniciales son simplemente:
    ${\vec {f}}^{1}={\vec {b}}^{1}={\begin{bmatrix}{1 \over M_{11}}\end{bmatrix}}={\begin{bmatrix}{1 \over t_{0}}\end{bmatrix}}$

### Usando los vectores de retroceso

Los pasos anteriores dan los N vectores de retroceso para M. De ahí, una ecuación más arbitraria es:
    ${\vec {y}}=\mathbf {M} \ {\vec {x}}$
La solución puede construirse de la misma manera iterativa en se construyeron los vectores de retroceso. Así, ${\vec {x}}$ debe generalizarse a la secuencia ${\vec {x}}^{n}$, de la cual ${\vec {x}}^{N}={\vec {x}}$.
La solución se construye entonces de manera iterativa, notando que si:
    $\mathbf {T} ^{n-1}{\begin{bmatrix}x_{1}^{n-1}\\x_{2}^{n-1}\\\dots \\x_{n-1}^{n-1}\\\end{bmatrix}}={\begin{bmatrix}y_{1}\\y_{2}\\\dots \\y_{n-1}\\\end{bmatrix}}$
Entonces, extendiendo con un cero nuevamente, y definiendo una constante de error donde corresponda:
    $\mathbf {T} ^{n}{\begin{bmatrix}x_{1}^{n-1}\\x_{2}^{n-1}\\\dots \\x_{n-1}^{n-1}\\0\\\end{bmatrix}}={\begin{bmatrix}y_{1}\\y_{2}\\\dots \\y_{n-1}\\\epsilon _{x}^{n-1}\end{bmatrix}}$
Podemos usar entonces el vector de retroceso n-ésimo para eliminar el término de error y reemplazarlo con la fórmula deseada de la siguiente manera:
    $\mathbf {T} ^{n}\left({\begin{bmatrix}x_{1}^{n-1}\\x_{2}^{n-1}\\\dots \\x_{n-1}^{n-1}\\0\\\end{bmatrix}}+(y_{n}-\epsilon _{x}^{n-1})\ {\vec {b}}^{n}\right)={\begin{bmatrix}y_{1}\\y_{2}\\\dots \\y_{n-1}\\y_{n}\\\end{bmatrix}}$
Extendiendo este método hasta n = N se obtiene la solución ${\vec {x}}$.
En la práctica, estos pasos se hacen a menudo en forma concurrente con el resto del procedimiento, pero forman una unidad coherente y merecen ser tratados como pasos separados.