Algoritmo de Levinson

De Wikipedia, la enciclopedia libre
Saltar a: navegación, búsqueda

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 Θ(n2), una mejora considerable frente a la eliminación de Gauss-Jordan, cuyo costo es de Θ(n3).

Hay algoritmos nuevos, llamados asintóticamente rápidos o a menudo algoritmos de Toeplitz superrápidos, que pueden resolver con un costo de Θ(n logpn) 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 [1]).

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 4n2 y luego a 3n2 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[editar]

Contexto[editar]

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 xi 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 Tn es una matriz de n×n que copia el bloque superior izquierdo de n×n de M, es decir, Tnij = Mij.

Tn 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[editar]

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 bni = fnn+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[editar]

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 Tn-1 a Tn, 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[editar]

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.

Referencias[editar]

  • Levinson, N. (1947). "The Wiener RMS error criterion in filter design and prediction," J. Math. Phys., v. 25, pp. 261-278.

Véase también[editar]