Recursión de Levinson

Algoritmo recursivo en álgebra lineal

La recursión de Levinson o recursión de Levinson-Durbin es un procedimiento de álgebra lineal para calcular recursivamente la solución de una ecuación que involucra una matriz de Toeplitz . El algoritmo se ejecuta en un tiempo Θ ( n ² ) , lo que es una gran mejora con respecto a la eliminación de Gauss-Jordan , que se ejecuta en Θ( n ³ ).

El algoritmo de Levinson-Durbin fue propuesto por primera vez por Norman Levinson en 1947, mejorado por James Durbin en 1960 y posteriormente mejorado a 4 n ² y luego a 3 n ² multiplicaciones por WF 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 con estas, la recursión de Levinson (en particular, la recursión de Levinson dividida) tiende a ser más rápida en términos computacionales, pero más sensible a imprecisiones computacionales como errores de redondeo .

El algoritmo de Bareiss para matrices de Toeplitz (que no debe confundirse con el algoritmo general de Bareiss ) se ejecuta aproximadamente tan rápido como la recursión de Levinson, pero utiliza el espacio O ( n ² ) , mientras que la recursión de Levinson utiliza solo el espacio O ( n ). Sin embargo, el algoritmo de Bareiss es numéricamente estable , ^{[1] [2]} mientras que la recursión de Levinson es, en el mejor de los casos, solo débilmente estable (es decir, exhibe estabilidad numérica para sistemas lineales bien condicionados ). ^[3]

Los algoritmos más nuevos, llamados algoritmos de Toeplitz asintóticamente rápidos o a veces superrápidos , pueden resolver en Θ( n log ^p n ) para varios p (por ejemplo, p = 2, ^{[4] [5]} p = 3 ^[6] ). La recursión de Levinson sigue siendo popular por varias razones; por un lado, es relativamente fácil de entender en comparación; por otro, puede ser más rápida que un algoritmo superrápido para n pequeños (generalmente n  < 256). ^[7]

Derivación

Fondo

Las ecuaciones matriciales siguen la forma

${\displaystyle \mathbf {M} \,{\vec {x}}={\vec {y}}.}$

El algoritmo de Levinson-Durbin puede utilizarse para cualquier ecuación de este tipo, siempre que M sea una matriz de Toeplitz conocida con una diagonal principal distinta de cero. Aquí se conoce un vector , y es un vector desconocido de números x _i aún por determinar. ${\displaystyle {\vec {y}}}$ ${\displaystyle {\vec {x}}}$

Para los fines de este artículo, ê _i es un vector formado enteramente por ceros, excepto por su i ésimo lugar, que contiene el valor uno. Su longitud estará determinada implícitamente por el contexto circundante. El término N se refiere al ancho de la matriz anterior: M es una matriz N × N. Finalmente, en este artículo, los superíndices se refieren a un índice inductivo , mientras que los subíndices denotan índices. Por ejemplo (y definición), en este artículo, la matriz T ⁿ es una matriz n × n que copia el bloque n × n superior izquierdo de M  , es decir, T ⁿ _ij = M _ij .

T ⁿ también es una matriz de Toeplitz, lo que significa que se puede escribir como

${\displaystyle \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 primer paso, se establecen dos conjuntos de vectores, llamados vectores directos y vectores inversos . Los vectores directos se utilizan para ayudar a obtener el conjunto de vectores inversos; luego, se pueden descartar inmediatamente. Los vectores inversos son necesarios para el segundo paso, donde se utilizan para construir la solución deseada.

La recursión de Levinson-Durbin define el n- ^ésimo "vector hacia adelante", denotado , como el vector de longitud n que satisface: ${\displaystyle {\vec {f}}^{n}}$

${\displaystyle \mathbf {T} ^{n}{\vec {f}}^{n}={\hat {e}}_{1}.}$

El n ^-ésimo "vector hacia atrás" se define de manera similar; es el vector de longitud n que satisface: ${\displaystyle {\vec {b}}^{n}}$

${\displaystyle \mathbf {T} ^{n}{\vec {b}}^{n}={\hat {e}}_{n}.}$

Una simplificación importante puede ocurrir cuando M es una matriz simétrica ; entonces los dos vectores están relacionados por b ⁿ _i = f ⁿ _{n +1− i} —es decir, son reversiones de fila entre sí. Esto puede ahorrar algunos cálculos adicionales en ese caso especial.

Obtención de los vectores hacia atrás

Incluso si la matriz no es simétrica, entonces el n ^-ésimo vector hacia adelante y hacia atrás se puede encontrar a partir de los vectores de longitud n  − 1 de la siguiente manera. Primero, el vector hacia adelante se puede extender con un cero para obtener:

${\displaystyle \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}}.}$

Al pasar de T ^{n −1} a T ⁿ , la columna adicional agregada a la matriz no perturba la solución cuando se utiliza un cero para extender el vector hacia delante. Sin embargo, la fila adicional agregada a la matriz ha perturbado la solución y ha creado un término de error no deseado ε _f que aparece en el último lugar. La ecuación anterior le da el valor de:

${\displaystyle \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}.}$

Volveremos a este error en breve y lo eliminaremos del nuevo vector hacia delante; pero primero, el vector hacia atrás debe extenderse de una manera similar (aunque invertida). Para el vector hacia atrás,

${\displaystyle \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}}.}$

Como antes, la columna adicional agregada a la matriz no altera este nuevo vector inverso, pero la fila adicional sí lo hace. Aquí tenemos otro error no deseado ε _b con valor:

${\displaystyle \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}.\ }$

Estos dos términos de error se pueden utilizar para formar vectores hacia adelante y hacia atrás de orden superior que se describen a continuación. Utilizando la linealidad de las matrices, la siguiente identidad se cumple para todos los : ${\displaystyle (\alpha ,\beta )}$

${\displaystyle \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}}.}$

Si se eligen α y β de modo que el lado derecho dé ê ₁ o ê _n , entonces la cantidad entre paréntesis cumplirá con la definición del n º vector ^hacia adelante o hacia atrás, respectivamente. Con esos alfa y beta elegidos, la suma de vectores entre paréntesis es simple y da el resultado deseado.

Para encontrar estos coeficientes, , son tales que: ${\displaystyle \alpha _{f}^{n}}$ ${\displaystyle \beta _{f}^{n}}$

${\displaystyle {\vec {f}}^{n}=\alpha _{f}^{n}{\begin{bmatrix}{\vec {f}}^{n-1}\\0\end{bmatrix}}+\beta _{f}^{n}{\begin{bmatrix}0\\{\vec {b}}^{n-1}\end{bmatrix}}}$

y respectivamente , son tales que: ${\displaystyle \alpha _{b}^{n}}$ ${\displaystyle \beta _{b}^{n}}$

${\displaystyle {\vec {b}}^{n}=\alpha _{b}^{n}{\begin{bmatrix}{\vec {f}}^{n-1}\\0\end{bmatrix}}+\beta _{b}^{n}{\begin{bmatrix}0\\{\vec {b}}^{n-1}\end{bmatrix}}.}$

Al multiplicar ambas ecuaciones anteriores por uno se obtiene la siguiente ecuación: ${\displaystyle {\mathbf {T} }^{n}}$

${\displaystyle {\begin{bmatrix}1&\varepsilon _{b}^{n}\\0&0\\\vdots &\vdots \\0&0\\\varepsilon _{f}^{n}&1\end{bmatrix}}{\begin{bmatrix}\alpha _{f}^{n}&\alpha _{b}^{n}\\\beta _{f}^{n}&\beta _{b}^{n}\end{bmatrix}}={\begin{bmatrix}1&0\\0&0\\\vdots &\vdots \\0&0\\0&1\end{bmatrix}}.}$

Ahora, descartando y colapsando todos los ceros en el medio de los dos vectores anteriores, solo queda la siguiente ecuación:

${\displaystyle {\begin{bmatrix}1&\varepsilon _{b}^{n}\\\varepsilon _{f}^{n}&1\end{bmatrix}}{\begin{bmatrix}\alpha _{f}^{n}&\alpha _{b}^{n}\\\beta _{f}^{n}&\beta _{b}^{n}\end{bmatrix}}={\begin{bmatrix}1&0\\0&1\end{bmatrix}}.}$

Una vez resueltos estos problemas (utilizando la fórmula inversa de la matriz 2×2 de Cramer), los nuevos vectores hacia adelante y hacia atrás son:

${\displaystyle {\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}}}$

${\displaystyle {\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}}.}$

Realizando estas sumas vectoriales, se obtienen los n ^-ésimos vectores hacia adelante y hacia atrás a partir de los anteriores. Todo lo que queda es encontrar el primero de estos vectores y, luego, algunas sumas y multiplicaciones rápidas dan los restantes. Los primeros vectores hacia adelante y hacia atrás son simplemente:

${\displaystyle {\vec {f}}^{1}={\vec {b}}^{1}=\left[{1 \over M_{11}}\right]=\left[{1 \over t_{0}}\right].}$

Usando los vectores hacia atrás

Los pasos anteriores dan los N vectores inversos para M. A partir de allí, una ecuación más arbitraria es:

${\displaystyle {\vec {y}}=\mathbf {M} \ {\vec {x}}.}$

La solución se puede construir de la misma manera recursiva en que se construyeron los vectores inversos. Por lo tanto, se debe generalizar a una secuencia de intermediarios , tales que . ${\displaystyle {\vec {x}}}$ ${\displaystyle {\vec {x}}^{n}}$ ${\displaystyle {\vec {x}}^{N}={\vec {x}}}$

Luego, la solución se construye de forma recursiva teniendo en cuenta que si

${\displaystyle \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}}.}$

Luego, extendiendo nuevamente con un cero y definiendo una constante de error donde sea necesario:

${\displaystyle \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}}.}$

Luego podemos utilizar el n ^-ésimo vector inverso para eliminar el término de error y reemplazarlo con la fórmula deseada de la siguiente manera:

${\displaystyle \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}}.}$

Extendiendo este método hasta que n = N obtenemos la solución . ${\displaystyle {\vec {x}}}$

En la práctica, estos pasos suelen realizarse simultáneamente con el resto del procedimiento, pero forman una unidad coherente y merecen ser tratados como un paso independiente.

Algoritmo de Levinson en bloque

Si M no es estrictamente Toeplitz, sino Toeplitz en bloque , la recursión de Levinson se puede derivar de la misma manera considerando la matriz de Toeplitz en bloque como una matriz de Toeplitz con elementos de matriz (Musicus 1988). Las matrices de Toeplitz en bloque surgen naturalmente en los algoritmos de procesamiento de señales cuando se trabaja con múltiples flujos de señales (por ejemplo, en sistemas MIMO ) o señales cicloestacionarias.

Véase también

• Recursión de Levinson dividida
• Predicción lineal
• Modelo autorregresivo

Notas

1. Bojanczyk y otros (1995).
2. Brent (1999).
3. Krishna y Wang (1993).
4. "Copia archivada" (PDF) . Archivado desde el original (PDF) el 25 de marzo de 2012. Consultado el 1 de abril de 2013 .{{cite web}}: CS1 maint: archived copy as title (link)
5. "Copia archivada" (PDF) . Archivado desde el original (PDF) el 15 de noviembre de 2009. Consultado el 28 de abril de 2009 .{{cite web}}: CS1 maint: archived copy as title (link)
6. "Copia archivada" (PDF) . saaz.cs.gsu.edu . Archivado desde el original (PDF) el 18 de abril de 2007 . Consultado el 12 de enero de 2022 .{{cite web}}: CS1 maint: archived copy as title (link)
7. "Copia archivada" (PDF) . Archivado desde el original (PDF) el 5 de septiembre de 2006. Consultado el 15 de agosto de 2006 .{{cite web}}: CS1 maint: archived copy as title (link)

Referencias

Definición de fuentes

• Levinson, N. (1947). "El criterio de error RMS de Wiener en el diseño y predicción de filtros". J. Math. Phys. , vol. 25, págs. 261–278.
• Durbin, J. (1960). "El ajuste de modelos de series temporales". Rev. Inst. Int. Stat. , v. 28, págs. 233–243.
• Trench, WF (1964). "Un algoritmo para la inversión de matrices de Toeplitz finitas". J. Soc. Indust. Appl. Math. , vol. 12, págs. 515–522.
• Musicus, BR (1988). "Algoritmos de Levinson y Fast Choleski para matrices de Toeplitz y casi Toeplitz". RLE TR No. 538, MIT. [1]
• Delsarte, P. y Genin, YV (1986). "El algoritmo Levinson dividido". IEEE Transactions on Acoustics, Speech, and Signal Processing , v. ASSP-34(3), págs. 470–478.

Trabajos futuros

• Bojanczyk, AW; Brent, RP; De Hoog, FR; Sweet, DR (1995). "Sobre la estabilidad de los algoritmos de factorización de Bareiss y Toeplitz relacionados". Revista SIAM sobre análisis de matrices y aplicaciones . 16 : 40–57. arXiv : 1004.5510 . doi :10.1137/S0895479891221563. S2CID  367586.
• Brent RP (1999), "Estabilidad de algoritmos rápidos para sistemas lineales estructurados", Algoritmos rápidos y confiables para matrices con estructura (editores—T. Kailath, AH Sayed), cap.4 ( SIAM ).
• Bunch, JR (1985). "Estabilidad de los métodos para resolver sistemas de ecuaciones de Toeplitz". SIAM J. Sci. Stat. Comput. , v. 6, págs. 349–364. [2]
• Krishna, H.; Wang, Y. (1993). "El algoritmo Split Levinson es débilmente estable". Revista SIAM sobre análisis numérico . 30 (5): 1498–1508. doi :10.1137/0730078.

Resúmenes

• Bäckström, T. (2004). "2.2. Recursión de Levinson–Durbin". Modelado predictivo lineal del habla: restricciones y descomposición de pares de líneas espectrales. Tesis doctoral. Informe n.º 71 / Universidad Tecnológica de Helsinki, Laboratorio de acústica y procesamiento de señales de audio. Espoo, Finlandia. [3]
• Claerbout, Jon F. (1976). "Capítulo 7 – Aplicaciones de formas de onda de los mínimos cuadrados". Fundamentos del procesamiento de datos geofísicos. Palo Alto: Blackwell Scientific Publications. [4]
• Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007), "Sección 2.8.2. Matrices de Toeplitz", Recetas numéricas: el arte de la computación científica (3.ª ed.), Nueva York: Cambridge University Press, ISBN 978-0-85-0-312-0 978-0-521-88068-8
• Golub, GH y Loan, CF Van (1996). "Sección 4.7: Toeplitz y sistemas relacionados" Matrix Computations , Johns Hopkins University Press