Algoritmo de Lanczos

El método iterativo para encontrar el mayor valor propio de una matriz ${\displaystyle A\,}$ puede resumirse señalando que si ${\displaystyle x_{0}\,}$ es un vector aleatorio y ${\displaystyle x_{n+1}=Ax_{n}\,}$ , entonces para un ${\displaystyle n}$ grande límite ${\displaystyle x_{n}/\|x_{n}\|}$ se acerca al vector propio normalizado correspondiente al valor propio más grande.

Si ${\displaystyle A=U\operatorname {diag} (\sigma _{i})U'\,}$ es la descomposición en valores propios de una matriz de ${\displaystyle A\,}$ , entonces ${\displaystyle A^{n}=U\operatorname {diag} (\sigma _{i}^{n})U'}$ . Como ${\displaystyle n\,}$ se hace grande, la matriz diagonal de valores propios estará acotada por el valor propio más grande(despreciando el caso en que el valor más grande este repetido ). Mientras, ${\displaystyle \|x_{n+1}\|/\|x_{n}\|\,}$ convergerán al mayor valor propio y ${\displaystyle x_{n}/\|x_{n}\|\,}$ al vector propio asociado. Si el mayor valor propio es múltiple, entonces ${\displaystyle x_{n}\,}$ convergerá a un vector en el subespacio generado por los vectores propios asociados con esos valores propios más grandes. Luego de haber encontrado el primer vector propio / valor, uno puede restringir sucesivamente el algoritmo para el espacio nulo (kernel) de los vectores propios conocidos para obtener el segundo mayor vector propio / valores y así sucesivamente.

En la práctica, este sencillo algoritmo no funciona muy bien para el cálculo de muchos de los vectores propios, ya que cualquier error de redondeo podría introducir ligeros cambios a los componentes de los vectores propios más significativos de nuevo en el cálculo, degradando la exactitud del cálculo.

Durante el procedimiento de aplicación del método, al obtener el último valor propio ${\displaystyle A^{n-1}v}$ , también se obtiene una serie de vectores ${\displaystyle A^{j}v,\,j=0,1,\cdots ,n-2}$ los cuales son descartados. Como ${\displaystyle n}$ es con frecuencia grande, puede que se tenga una gran cantidad de información que no se utilice. Los algoritmos más avanzados, Como el el algoritmo de Arnoldi, el algoritmo de Lanczos, guardan esta información y utilizan el proceso de Gram–Schmidt o el algoritmo Householder para ortogonalizar nuevamente en una base que abarca el subespacio de Krylov correspondiente a la matriz ${\displaystyle A}$ .

El algoritmo

El algoritmo de Lanczos puede ser visto como un sistema simplificado del algoritmo de Arnoldi que se aplica a matrices hermitianas. El ${\displaystyle m-esimo}$ paso del algoritmo transforma la matriz ${\displaystyle A}$ en una matriz tridiagonal ${\displaystyle T_{mm}}$ ; cuando ${\displaystyle m}$ es igual a la dimensión de ${\displaystyle A}$ , ${\displaystyle T_{mm}}$ es semejante a ${\displaystyle A}$ .

Definiciones

Se quiere calcular la matriz tridiagonal y simétrica ${\displaystyle T_{mm}=V_{m}^{*}AV_{m}.}$

Los elementos de la diagonal se denotan por ${\displaystyle \alpha _{j}=t_{jj}}$ , Y los elementos fuera de la diagonal son denotados por ${\displaystyle \beta _{j}=t_{j-1,j}}$ .

Notar que ${\displaystyle t_{j-1,j}=t_{j,j-1}}$ , debido a la simetría de la matriz.

Iteración

(Nota: siguiendo estos pasos solamente no se encontraran correctamente los valores y vectores propios. Se deben tener en cuenta más consideraciones para corregir errores numéricos. Ver la sección Estabilidad numérica.)

En principio, hay cuatro maneras de escribir el procedimiento de iteración. Paige[1972] y otros trabajos muestran que el siguiente procedimiento es el más estable numéricamente.^[6]​^[7]​

 ${\displaystyle v_{1}\leftarrow \,}$  vector aleatorio con norma 1. ${\displaystyle v_{0}\leftarrow 0\,}$  ${\displaystyle \beta _{1}\leftarrow 0\,}$ 

 for ${\displaystyle j=1,2,\cdots ,m-1\,}$    ${\displaystyle w_{j}\leftarrow Av_{j}\,}$    ${\displaystyle \alpha _{j}\leftarrow w_{j}\cdot v_{j}\,}$    ${\displaystyle w_{j}\leftarrow w_{j}-\alpha _{j}v_{j}-\beta _{j}v_{j-1}\,}$    ${\displaystyle \beta _{j+1}\leftarrow \left\|w_{j}\right\|\,}$    ${\displaystyle v_{j+1}\leftarrow w_{j}/\beta _{j+1}\,}$  endfor

 ${\displaystyle w_{m}\leftarrow Av_{m}\,}$  ${\displaystyle \alpha _{m}\leftarrow w_{m}\cdot v_{m}\,}$  return

Aquí, ${\displaystyle x\cdot y}$ representa el producto escalar de vectores ${\displaystyle x}$ y ${\displaystyle y}$ .

Después de la iteración, se obtiene la ${\displaystyle \alpha _{j}}$ and ${\displaystyle \beta _{j}}$ con los que se forma una matriz tridiagonal.

${\displaystyle T_{mm}={\begin{pmatrix}\alpha _{1}&\beta _{2}&&&&0\\\beta _{2}&\alpha _{2}&\beta _{3}&&&\\&\beta _{3}&\alpha _{3}&\ddots &&\\&&\ddots &\ddots &\beta _{m-1}&\\&&&\beta _{m-1}&\alpha _{m-1}&\beta _{m}\\0&&&&\beta _{m}&\alpha _{m}\\\end{pmatrix}}}$

Los vectores ${\displaystyle v_{j}}$ (vectores de Lanczos) generados sobre la marcha forman la matriz de transformación${\displaystyle V_{m}=\left(v_{1},v_{2},\cdots ,v_{m}\right)}$ , que es útil para el cálculo de los vectores propios (ver más abajo). En la práctica, se podría guardar en cada iteración (pero ocupa memoria), o se podría recalcular cuando se necesite, siempre y cuando se tenga el primer vector ${\displaystyle v_{1}}$ . En cada iteración del algoritmo realiza una multiplicación matriz-vector y otras operaciones de punto flotantes.

Encontrando valores y vectores propios

Después de que la matriz ${\displaystyle T_{mm}}$ es calculada, se pueden encontrar sus valores propios ${\displaystyle \lambda _{i}^{(m)}}$ y sus correspondientes vectores propios ${\displaystyle u_{i}^{(m)}}$ (Por ejemplo, usando el algoritmo QR o Multiple Relatively Robust Representations (MRRR)). Los valores y vectores propios de${\displaystyle T}$ se pueden obtener en ${\displaystyle {\mathcal {O}}(m^{2})}$ usando MRRR; si se quiere buscar solo los valores propios estos se calculan en ${\displaystyle {\mathcal {O}}(m^{2})}$ usando bisección espectral.

Se puede demostrar que los valores propios son valores propios aproximados de la matriz original ${\displaystyle A}$ .

Los vectores propios ${\displaystyle y_{i}}$ de ${\displaystyle A}$ se pueden calcular ${\displaystyle y_{i}=V_{m}u_{i}^{(m)}}$ , donde ${\displaystyle V_{m}}$ es la matriz de transformación cuyos vectores columna son ${\displaystyle v_{1},v_{2},\cdots ,v_{m}}$ .

Estabilidad Numérica

Estabilidad significa cuánto se verá afectado el algoritmo (es decir, se encontrara un resultado aproximado cercano al original) si hay pequeños errores numéricos introducidos y acumulados. Estabilidad numérica es el criterio central para juzgar la utilidad de una implementación de un algoritmo en un ordenador con redondeo.Para el algoritmo de Lanczos, se puede demostrar que con una aritmética exacta, el conjunto de vectores ${\displaystyle v_{1},v_{2},\cdots ,v_{m+1}}$ forman una base ortonormal, los valores y vectores propios encontrados son buenas aproximaciones a los de la matriz original. Sin embargo, en la práctica (como los cálculos se realizan en aritméticas de punto flotante donde la inexactitud es inevitable, la ortogonalidad se pierde rápidamente y en algunos casos el nuevo vector incluso podría ser linealmente dependiente del conjunto que ya está construido. Como resultado, algunos de los valores propios de la matriz tridiagonal resultante no son buenas aproximaciones para los de la matriz original. Por lo que, el algoritmo de Lanczos no es muy estable.

Los que usen el algoritmo deben ser capaz de encontrar y eliminar los valores propios "con errores". Las implementaciones prácticas del algoritmo de Lanczos van en tres direcciones para luchar contra este problema de estabilidad:^[6]​^[7]​

1. Prevenir la pérdida de ortogonalidad
2. Recuperar la ortogonalidad después de que se genera la base.
3. Después de identificar los valores propios "con errores", eliminarlos.

Existen variantes en el algoritmo de Lanczos donde los vectores involucrados son vectores columnas, matrices estrechas y las constantes de la normalización son pequeñas matrices cuadradas. Estos se llaman "bloques" y el algoritmo de Lanczos puede ser mucho más rápido en computadoras con un gran número de registros y tiempos de esfuerzo de memoria largos.

Muchas implementaciones del algoritmo de Lanczos reinician después de un cierto número de iteraciones. Una de las variaciones del reinicio más influyentes es el método de Lanczos con reinicio implícito,^[8]​ que se implementa enARPACK.^[9]​ Esto ha llevado a un número de otras variantes de reinicio tales como reinicio bidiagonalizado de Lanczos.^[10]​ Otra variante de reinicio exitosa es el método Thick-Restart de Lanczos,^[11]​ que se ha implementado en un paquete de software llamado TRLAN.^[12]​

Kernel sobre un cuerpo finito

En 1995, Peter Montgomery publicó un algoritmo, basado en el algoritmo de Lanczos para encontrar elementos del kernel de una gran matriz sparse sobre GF(2); ya que el conjunto de personas interesadas en matrices sparse de gran dimensión sobre cuerpos finitos y el conjunto de personas interesadas en los problemas de los valores propios más grandes, tienen pocas personas en común, esto es a menudo también llamado el algoritmo de Lanczos por bloques.

Los algoritmos Lanczos son muy atractivos porque la multiplicación por ${\displaystyle A\,}$ es la única operación lineal a gran escala. Dado que los motores de recuperación de texto implementan esta operación, el algoritmo Lanczos se puede aplicar de manera eficiente a los documentos de texto (ver indexación semántica latente). Los vectores propios son también importantes para los métodos de clasificación a gran escala tales como el algoritmo HITS desarrollado por Jon Kleinberg, o el PageRank algoritmo usado por Google.

Algoritmos de Lanczos también se utilizan en Física de la materia condensada como un método para resolver hamiltonianas de sistemas de electrones fuertemente correlacionados.^[13]​

La biblioteca NAG contiene varias métodos^[14]​ para la solución de sistemas lineales grandes y problemas de valores/vectores propios que utilizan el algoritmo de Lanczos.

MATLAB y GNU Octave vienen con ARPACK incorporado. Ambos analizan matrices usando la función eigs() (Matlab/Octave).

La implementación de Matlab del algoritmo de Lanczos (nota: problemas de precisión) está disponible como parte de la Gaussian Belief Propagation Matlab Package. El GraphLab^[15]​ biblioteca de filtrado colaborativo incorpora una implementación paralelizada del algoritmo de Lanczos (en C ++) para multinúcleo.La biblioteca PRIMME también implementa el algoritmo de Lanczos.

• Golub and van Loan give very good descriptions of the various forms of Lanczos algorithms in their book Matrix Computations
• Andrew Ng et al., an analysis of PageRank
• Lanczos and conjugate gradient methods B. A. LaMacchia and A. M. Odlyzko, Solving Large Sparse Linear Systems Over Finite Fields.