Algoritmo de Lanczos

Cálculo de valores propios numéricos

El algoritmo de Lanczos es un método iterativo ideado por Cornelius Lanczos que es una adaptación de los métodos de potencia para encontrar los valores propios y vectores propios "más útiles" (que tienden hacia los extremos más altos/más bajos) de una matriz hermítica , donde es a menudo, pero no necesariamente, mucho menor que . ^[1] Aunque computacionalmente eficiente en principio, el método tal como se formuló inicialmente no era útil, debido a su inestabilidad numérica . ${\displaystyle m}$ ${\displaystyle n\times n}$ ${\displaystyle m}$ ${\displaystyle n}$

En 1970, Ojalvo y Newman demostraron cómo hacer que el método fuera numéricamente estable y lo aplicaron a la solución de estructuras de ingeniería muy grandes sometidas a carga dinámica. ^[2] Esto se logró utilizando un método para purificar los vectores de Lanczos (es decir, reortogonalizando repetidamente cada vector recién generado con todos los generados previamente) ^[2] hasta cualquier grado de precisión, lo que cuando no se realizó, produjo una serie de vectores que estaban altamente contaminados por aquellos asociados con las frecuencias naturales más bajas.

En su trabajo original, estos autores también sugirieron cómo seleccionar un vector de inicio (es decir, usar un generador de números aleatorios para seleccionar cada elemento del vector de inicio) y sugirieron un método determinado empíricamente para determinar , el número reducido de vectores (es decir, debe seleccionarse para que sea aproximadamente 1,5 veces el número de valores propios precisos deseados). Poco después, su trabajo fue seguido por Paige, quien también proporcionó un análisis de errores. ^{[3] [4]} En 1988, Ojalvo produjo una historia más detallada de este algoritmo y una prueba de error de valor propio eficiente. ^[5] ${\displaystyle m}$

El algoritmo

Introduzca una matriz hermítica de tamaño y, opcionalmente, un número de iteraciones (por defecto, sea ). ${\displaystyle A}$ ${\displaystyle n\times n}$ ${\displaystyle m}$ ${\displaystyle m=n}$

• En sentido estricto, el algoritmo no necesita acceder a la matriz explícita, sino únicamente a una función que calcule el producto de la matriz por un vector arbitrario. Esta función se invoca en la mayoría de los casos. ${\displaystyle v\mapsto Av}$ ${\displaystyle m}$

Genere una matriz con columnas ortonormales y una matriz simétrica real tridiagonal de tamaño . Si , entonces es unitaria y . ${\displaystyle n\times m}$ ${\displaystyle V}$ ${\displaystyle T=V^{*}AV}$ ${\displaystyle m\times m}$ ${\displaystyle m=n}$ ${\displaystyle V}$ ${\displaystyle A=VTV^{*}}$

Advertencia La iteración de Lanczos es propensa a la inestabilidad numérica. Cuando se ejecuta en aritmética no exacta, se deben tomar medidas adicionales (como se describe en secciones posteriores) para garantizar la validez de los resultados.

1. Sea un vector arbitrario con norma euclidiana . ${\displaystyle v_{1}\in \mathbb {C} ^{n}}$ ${\displaystyle 1}$
2. Paso de iteración inicial abreviado:
   1. Dejar . ${\displaystyle w_{1}'=Av_{1}}$
   2. Dejar . ${\displaystyle \alpha _{1}=w_{1}'^{*}v_{1}}$
   3. Dejar . ${\displaystyle w_{1}=w_{1}'-\alpha _{1}v_{1}}$
3. Para hacer: ${\displaystyle j=2,\dots ,m}$
   1. Sea (también norma euclidiana ). ${\displaystyle \beta _{j}=\|w_{j-1}\|}$
   2. Si , entonces sea , ${\displaystyle \beta _{j}\neq 0}$ ${\displaystyle v_{j}=w_{j-1}/\beta _{j}}$

      de lo contrario, elija un vector arbitrario con norma euclidiana que sea ortogonal a todos los . ${\displaystyle v_{j}}$ ${\displaystyle 1}$ ${\displaystyle v_{1},\dots ,v_{j-1}}$

   3. Dejar . ${\displaystyle w_{j}'=Av_{j}}$
   4. Dejar . ${\displaystyle \alpha _{j}=w_{j}'^{*}v_{j}}$
   5. Dejar . ${\displaystyle w_{j}=w_{j}'-\alpha _{j}v_{j}-\beta _{j}v_{j-1}}$
4. Sea la matriz con columnas . Sea . ${\displaystyle V}$ ${\displaystyle v_{1},\dots ,v_{m}}$ ${\displaystyle T={\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}}}$

Nota para . ${\displaystyle Av_{j}=w_{j}'=\beta _{j+1}v_{j+1}+\alpha _{j}v_{j}+\beta _{j}v_{j-1}}$ ${\displaystyle 2<j<m}$

En principio, existen cuatro formas de escribir el procedimiento de iteración. Paige y otros trabajos muestran que el orden de operaciones anterior es el más estable numéricamente. ^{[6] [7]} En la práctica, el vector inicial puede tomarse como otro argumento del procedimiento, y se pueden incluir indicadores de imprecisión numérica como condiciones adicionales de terminación del bucle. ${\displaystyle v_{1}}$ ${\displaystyle \beta _{j}=0}$

Sin contar la multiplicación de matriz-vector, cada iteración realiza operaciones aritméticas. La multiplicación de matriz-vector se puede realizar en operaciones aritméticas donde es el número promedio de elementos distintos de cero en una fila. La complejidad total es, por lo tanto , o si ; el algoritmo de Lanczos puede ser muy rápido para matrices dispersas. Los esquemas para mejorar la estabilidad numérica generalmente se evalúan en función de este alto rendimiento. ${\displaystyle O(n)}$ ${\displaystyle O(dn)}$ ${\displaystyle d}$ ${\displaystyle O(dmn)}$ ${\displaystyle O(dn^{2})}$ ${\displaystyle m=n}$

Los vectores se denominan vectores de Lanczos . El vector no se utiliza después de que se calcula y el vector no se utiliza después de que se calcula. Por lo tanto, se puede utilizar el mismo almacenamiento para los tres. Del mismo modo, si solo se busca la matriz tridiagonal, entonces la iteración en bruto no necesita después de haber calculado , aunque algunos esquemas para mejorar la estabilidad numérica lo necesitarían más adelante. A veces, los vectores de Lanczos posteriores se vuelven a calcular cuando es necesario. ${\displaystyle v_{j}}$ ${\displaystyle w_{j}'}$ ${\displaystyle w_{j}}$ ${\displaystyle w_{j}}$ ${\displaystyle v_{j+1}}$ ${\displaystyle T}$ ${\displaystyle v_{j-1}}$ ${\displaystyle w_{j}}$ ${\displaystyle v_{1}}$

Aplicación al problema propio

El algoritmo de Lanczos se utiliza con mayor frecuencia en el contexto de la búsqueda de los valores y vectores propios de una matriz, pero mientras que una diagonalización ordinaria de una matriz haría que los vectores y valores propios fueran evidentes a partir de la inspección, no sucede lo mismo con la tridiagonalización realizada por el algoritmo de Lanczos; se necesitan pasos adicionales no triviales para calcular incluso un único valor o vector propio. No obstante, la aplicación del algoritmo de Lanczos suele suponer un avance significativo en el cálculo de la descomposición propia.

Si es un valor propio de , y su vector propio ( ), entonces es un vector propio correspondiente de con el mismo valor propio: ${\displaystyle \lambda }$ ${\displaystyle T}$ ${\displaystyle x}$ ${\displaystyle Tx=\lambda x}$ ${\displaystyle y=Vx}$ ${\displaystyle A}$

$${\displaystyle {\begin{aligned}Ay&=AVx\\&=VTV^{*}Vx\\&=VTIx\\&=VTx\\&=V(\lambda x)\\&=\lambda Vx\\&=\lambda y.\end{aligned}}}$$

De esta forma, el algoritmo de Lanczos transforma el problema de descomposición propia para en el problema de descomposición propia para . ${\displaystyle A}$ ${\displaystyle T}$

1. Para matrices tridiagonales, existen varios algoritmos especializados, a menudo con una complejidad computacional mejor que los algoritmos de propósito general. Por ejemplo, si es una matriz simétrica tridiagonal, entonces: ${\displaystyle T}$ ${\displaystyle m\times m}$
   • La recursión continua permite calcular el polinomio característico en las operaciones y evaluarlo en un punto de las operaciones. ${\displaystyle O(m^{2})}$ ${\displaystyle O(m)}$
   • El algoritmo de valores propios de dividir y vencer se puede utilizar para calcular toda la descomposición propia de las operaciones. ${\displaystyle T}$ ${\displaystyle O(m^{2})}$
   • El método multipolar rápido ^[8] puede calcular todos los valores propios en sólo operaciones. ${\displaystyle O(m\log m)}$
2. Se sabe que algunos algoritmos generales de descomposición de vectores propios, en particular el algoritmo QR , convergen más rápido para matrices tridiagonales que para matrices generales. La complejidad asintótica del algoritmo QR tridiagonal es igual a la del algoritmo divide y vencerás (aunque el factor constante puede ser diferente); dado que los vectores propios juntos tienen elementos, esto es asintóticamente óptimo. ${\displaystyle O(m^{2})}$ ${\displaystyle m^{2}}$
3. Incluso los algoritmos cuyas tasas de convergencia no se ven afectadas por las transformaciones unitarias, como el método de potencia y la iteración inversa , pueden disfrutar de beneficios de rendimiento de bajo nivel al ser aplicados a la matriz tridiagonal en lugar de a la matriz original . Dado que es muy esparcido con todos los elementos distintos de cero en posiciones altamente predecibles, permite un almacenamiento compacto con un rendimiento excelente en comparación con el almacenamiento en caché . Del mismo modo, es una matriz real con todos los vectores propios y valores propios reales, mientras que en general puede tener elementos complejos y vectores propios, por lo que la aritmética real es suficiente para encontrar los vectores propios y valores propios de . ${\displaystyle T}$ ${\displaystyle A}$ ${\displaystyle T}$ ${\displaystyle T}$ ${\displaystyle A}$ ${\displaystyle T}$
4. Si es muy grande, entonces reducirlo para que tenga un tamaño manejable aún permitirá encontrar los valores propios y vectores propios más extremos de ; en la región, el algoritmo de Lanczos puede verse como un esquema de compresión con pérdida para matrices hermíticas, que enfatiza la preservación de los valores propios extremos. ${\displaystyle n}$ ${\displaystyle m}$ ${\displaystyle T}$ ${\displaystyle A}$ ${\displaystyle m\ll n}$

La combinación de un buen rendimiento para matrices dispersas y la capacidad de calcular varios valores propios (sin calcular todos) son las principales razones para elegir utilizar el algoritmo de Lanczos.

Aplicación a la tridiagonalización

Aunque el problema propio es a menudo la motivación para aplicar el algoritmo de Lanczos, la operación que el algoritmo realiza principalmente es la tridiagonalización de una matriz, para la cual se han favorecido las transformaciones de Householder numéricamente estables desde la década de 1950. Durante la década de 1960, el algoritmo de Lanczos fue ignorado. El interés en él fue rejuvenecido por la teoría de convergencia de Kaniel-Paige y el desarrollo de métodos para prevenir la inestabilidad numérica, pero el algoritmo de Lanczos sigue siendo el algoritmo alternativo que se prueba solo si Householder no es satisfactorio. ^[9]

Los aspectos en los que difieren ambos algoritmos incluyen:

• Lanczos aprovecha el hecho de ser una matriz dispersa, mientras que Householder no, y generará un relleno . ${\displaystyle A}$
• Lanczos trabaja todo el tiempo con la matriz original (y no tiene problemas con que ésta sea conocida sólo implícitamente), mientras que Householder en su versión básica quiere modificar la matriz durante el cálculo (aunque eso se puede evitar). ${\displaystyle A}$
• Cada iteración del algoritmo de Lanczos produce otra columna de la matriz de transformación final , mientras que una iteración de Householder produce otro factor en una factorización unitaria de . Sin embargo, cada factor está determinado por un único vector, por lo que los requisitos de almacenamiento son los mismos para ambos algoritmos y se pueden calcular a tiempo. ${\displaystyle V}$ ${\displaystyle Q_{1}Q_{2}\dots Q_{n}}$ ${\displaystyle V}$ ${\displaystyle V=Q_{1}Q_{2}\dots Q_{n}}$ ${\displaystyle O(n^{3})}$
• El jefe de familia es numéricamente estable, mientras que el Lanczos en bruto no lo es.
• Lanczos es altamente paralelo, con solo puntos de sincronización (los cálculos de y ). Householder es menos paralelo, ya que tiene una secuencia de cantidades escalares calculadas que dependen cada una de la cantidad anterior en la secuencia. ${\displaystyle O(n)}$ ${\displaystyle \alpha _{j}}$ ${\displaystyle \beta _{j}}$ ${\displaystyle O(n^{2})}$

Derivación del algoritmo

Hay varias líneas de razonamiento que conducen al algoritmo de Lanczos.

Un método de poder más previsor

El método de potencia para encontrar el valor propio de mayor magnitud y un vector propio correspondiente de una matriz es aproximadamente ${\displaystyle A}$

1. Elija un vector aleatorio . ${\displaystyle u_{1}\neq 0}$
2. Para (hasta que la dirección de haya convergido) hacer: ${\displaystyle j\geqslant 1}$ ${\displaystyle u_{j}}$
   1. Dejar ${\displaystyle u_{j+1}'=Au_{j}.}$
   2. Dejar ${\displaystyle u_{j+1}=u_{j+1}'/\|u_{j+1}'\|.}$

• En el límite grande, se aproxima al vector propio normado correspondiente al valor propio de mayor magnitud. ${\displaystyle j}$ ${\displaystyle u_{j}}$

Una crítica que se puede hacer a este método es que es un desperdicio: gasta mucho trabajo (los productos matriz-vector en el paso 2.1) extrayendo información de la matriz , pero presta atención solo al último resultado; las implementaciones suelen utilizar la misma variable para todos los vectores , haciendo que cada nueva iteración sobrescriba los resultados de la anterior. Puede ser conveniente, en cambio, conservar todos los resultados intermedios y organizar los datos. ${\displaystyle A}$ ${\displaystyle u_{j}}$

Una pieza de información que está disponible de manera trivial a partir de los vectores es una cadena de subespacios de Krylov . Una forma de afirmar esto sin introducir conjuntos en el algoritmo es afirmar que calcula ${\displaystyle u_{j}}$

un subconjunto de una base tal que para cada uno y todos ${\displaystyle \{v_{j}\}_{j=1}^{m}}$ ${\displaystyle \mathbb {C} ^{n}}$ ${\displaystyle Ax\in \operatorname {span} (v_{1},\dotsc ,v_{j+1})}$ ${\displaystyle x\in \operatorname {span} (v_{1},\dotsc ,v_{j})}$ ${\displaystyle 1\leqslant j<m;}$

esto se satisface trivialmente con siempre que sea linealmente independiente de (y en el caso de que exista tal dependencia, entonces se puede continuar la secuencia eligiendo como un vector arbitrario linealmente independiente de ). Sin embargo, es probable que una base que contenga los vectores esté numéricamente mal condicionada , ya que esta secuencia de vectores está diseñada para converger a un vector propio de . Para evitarlo, se puede combinar la iteración de potencia con un proceso de Gram-Schmidt , para producir en su lugar una base ortonormal de estos subespacios de Krylov. ${\displaystyle v_{j}=u_{j}}$ ${\displaystyle u_{j}}$ ${\displaystyle u_{1},\dotsc ,u_{j-1}}$ ${\displaystyle v_{j}}$ ${\displaystyle u_{1},\dotsc ,u_{j-1}}$ ${\displaystyle u_{j}}$ ${\displaystyle A}$

1. Elija un vector aleatorio de norma euclidiana . Sea . ${\displaystyle u_{1}}$ ${\displaystyle 1}$ ${\displaystyle v_{1}=u_{1}}$
2. Para hacer: ${\displaystyle j=1,\dotsc ,m-1}$
   1. Dejar . ${\displaystyle u_{j+1}'=Au_{j}}$
   2. Para todo sea . (Estas son las coordenadas de con respecto a los vectores base .) ${\displaystyle k=1,\dotsc ,j}$ ${\displaystyle g_{k,j}=v_{k}^{*}u_{j+1}'}$ ${\displaystyle Au_{j}=u_{j+1}'}$ ${\displaystyle v_{1},\dotsc ,v_{j}}$
   3. Sea . (Cancela el componente de que está en .) ${\displaystyle w_{j+1}=u_{j+1}'-\sum _{k=1}^{j}g_{k,j}v_{k}}$ ${\displaystyle u_{j+1}'}$ ${\displaystyle \operatorname {span} (v_{1},\dotsc ,v_{j})}$
   4. Si entonces sea y , ${\displaystyle w_{j+1}\neq 0}$ ${\displaystyle u_{j+1}=u_{j+1}'/\|u_{j+1}'\|}$ ${\displaystyle v_{j+1}=w_{j+1}/\|w_{j+1}\|}$

      de lo contrario, elija como un vector arbitrario de norma euclidiana que sea ortogonal a todos los . ${\displaystyle u_{j+1}=v_{j+1}}$ ${\displaystyle 1}$ ${\displaystyle v_{1},\dotsc ,v_{j}}$

La relación entre los vectores de iteración de potencia y los vectores ortogonales es que ${\displaystyle u_{j}}$ ${\displaystyle v_{j}}$

${\displaystyle Au_{j}=\|u_{j+1}'\|u_{j+1}=u_{j+1}'=w_{j+1}+\sum _{k=1}^{j}g_{k,j}v_{k}=\|w_{j+1}\|v_{j+1}+\sum _{k=1}^{j}g_{k,j}v_{k}}$.

Aquí se puede observar que en realidad no necesitamos los vectores para calcular estos , porque y por lo tanto la diferencia entre y está en , que se cancela mediante el proceso de ortogonalización. Por lo tanto, la misma base para la cadena de subespacios de Krylov se calcula mediante ${\displaystyle u_{j}}$ ${\displaystyle v_{j}}$ ${\displaystyle u_{j}-v_{j}\in \operatorname {span} (v_{1},\dotsc ,v_{j-1})}$ ${\displaystyle u_{j+1}'=Au_{j}}$ ${\displaystyle w_{j+1}'=Av_{j}}$ ${\displaystyle \operatorname {span} (v_{1},\dotsc ,v_{j})}$

1. Elija un vector aleatorio de norma euclidiana . ${\displaystyle v_{1}}$ ${\displaystyle 1}$
2. Para hacer: ${\displaystyle j=1,\dotsc ,m-1}$
   1. Dejar . ${\displaystyle w_{j+1}'=Av_{j}}$
   2. Para todos dejar . ${\displaystyle k=1,\dotsc ,j}$ ${\displaystyle h_{k,j}=v_{k}^{*}w_{j+1}'}$
   3. Dejar . ${\displaystyle w_{j+1}=w_{j+1}'-\sum _{k=1}^{j}h_{k,j}v_{k}}$
   4. Dejar . ${\displaystyle h_{j+1,j}=\|w_{j+1}\|}$
   5. Si entonces dejemos , ${\displaystyle h_{j+1,j}\neq 0}$ ${\displaystyle v_{j+1}=w_{j+1}/h_{j+1,j}}$

      de lo contrario, elija como un vector arbitrario de norma euclidiana que sea ortogonal a todos los . ${\displaystyle v_{j+1}}$ ${\displaystyle 1}$ ${\displaystyle v_{1},\dotsc ,v_{j}}$

A priori los coeficientes satisfacen ${\displaystyle h_{k,j}}$

${\displaystyle Av_{j}=\sum _{k=1}^{j+1}h_{k,j}v_{k}}$para todos ; ${\displaystyle j<m}$

La definición puede parecer un poco extraña, pero se ajusta al patrón general ya que ${\displaystyle h_{j+1,j}=\|w_{j+1}\|}$ ${\displaystyle h_{k,j}=v_{k}^{*}w_{j+1}'}$

${\displaystyle v_{j+1}^{*}w_{j+1}'=v_{j+1}^{*}w_{j+1}=\|w_{j+1}\|v_{j+1}^{*}v_{j+1}=\|w_{j+1}\|.}$

Debido a que los vectores de iteración de potencia que se eliminaron de esta recursión satisfacen los vectores y los coeficientes que contienen suficiente información para que se puedan calcular todos , no se perdió nada al cambiar los vectores. (De hecho, resulta que los datos recopilados aquí brindan aproximaciones significativamente mejores del valor propio más grande que las que se obtienen a partir de un número igual de iteraciones en el método de potencia, aunque eso no es necesariamente obvio en este punto). ${\displaystyle u_{j}}$ ${\displaystyle u_{j}\in \operatorname {span} (v_{1},\ldots ,v_{j}),}$ ${\displaystyle \{v_{j}\}_{j=1}^{m}}$ ${\displaystyle h_{k,j}}$ ${\displaystyle A}$ ${\displaystyle u_{1},\ldots ,u_{m}}$

Este último procedimiento es la iteración de Arnoldi . El algoritmo de Lanczos surge entonces como la simplificación que se obtiene al eliminar los pasos de cálculo que resultan triviales cuando es hermítico, en particular la mayoría de los coeficientes resultan ser cero. ${\displaystyle A}$ ${\displaystyle h_{k,j}}$

Elementalmente, si es hermitiano entonces ${\displaystyle A}$

${\displaystyle h_{k,j}=v_{k}^{*}w_{j+1}'=v_{k}^{*}Av_{j}=v_{k}^{*}A^{*}v_{j}=(Av_{k})^{*}v_{j}.}$

Porque sabemos que , y dado que por construcción es ortogonal a este subespacio, este producto interno debe ser cero. (Esta es esencialmente también la razón por la que a las sucesiones de polinomios ortogonales siempre se les puede dar una relación de recurrencia de tres términos ). Porque se obtiene ${\displaystyle k<j-1}$ ${\displaystyle Av_{k}\in \operatorname {span} (v_{1},\ldots ,v_{j-1})}$ ${\displaystyle v_{j}}$ ${\displaystyle k=j-1}$

${\displaystyle h_{j-1,j}=(Av_{j-1})^{*}v_{j}={\overline {v_{j}^{*}Av_{j-1}}}={\overline {h_{j,j-1}}}=h_{j,j-1}}$

ya que este último es real por ser la norma de un vector. Porque se obtiene ${\displaystyle k=j}$

${\displaystyle h_{j,j}=(Av_{j})^{*}v_{j}={\overline {v_{j}^{*}Av_{j}}}={\overline {h_{j,j}}},}$

lo que significa que esto también es real.

De manera más abstracta, si es la matriz con columnas , entonces los números pueden identificarse como elementos de la matriz , y para la matriz es Hessenberg superior . Dado que ${\displaystyle V}$ ${\displaystyle v_{1},\ldots ,v_{m}}$ ${\displaystyle h_{k,j}}$ ${\displaystyle H=V^{*}AV}$ ${\displaystyle h_{k,j}=0}$ ${\displaystyle k>j+1;}$ ${\displaystyle H}$

${\displaystyle H^{*}=\left(V^{*}AV\right)^{*}=V^{*}A^{*}V=V^{*}AV=H}$

La matriz es hermítica. Esto implica que también es Hessenberg inferior, por lo que, de hecho, debe ser tridiagional. Al ser hermítica, su diagonal principal es real y, dado que su primera subdiagonal es real por construcción, lo mismo es cierto para su primera superdiagonal. Por lo tanto, es una matriz real y simétrica: la matriz de la especificación del algoritmo de Lanczos. ${\displaystyle H}$ ${\displaystyle H}$ ${\displaystyle H}$ ${\displaystyle T}$

Aproximación simultánea de valores propios extremos

Una forma de caracterizar los vectores propios de una matriz hermítica es como puntos estacionarios del cociente de Rayleigh. ${\displaystyle A}$

${\displaystyle r(x)={\frac {x^{*}Ax}{x^{*}x}},\qquad x\in \mathbb {C} ^{n}.}$

En particular, el valor propio más grande es el máximo global de y el valor propio más pequeño es el mínimo global de . ${\displaystyle \lambda _{\max }}$ ${\displaystyle r}$ ${\displaystyle \lambda _{\min }}$ ${\displaystyle r}$

Dentro de un subespacio de baja dimensión de puede ser factible localizar el máximo y el mínimo de . Repitiendo eso para una cadena creciente se producen dos secuencias de vectores: y tales que y ${\displaystyle {\mathcal {L}}}$ ${\displaystyle \mathbb {C} ^{n}}$ ${\displaystyle x}$ ${\displaystyle y}$ ${\displaystyle r}$ ${\displaystyle {\mathcal {L}}_{1}\subset {\mathcal {L}}_{2}\subset \cdots }$ ${\displaystyle x_{1},x_{2},\ldots }$ ${\displaystyle y_{1},y_{2},\dotsc }$ ${\displaystyle x_{j},y_{j}\in {\mathcal {L}}_{j}}$

${\displaystyle {\begin{aligned}r(x_{1})&\leqslant r(x_{2})\leqslant \cdots \leqslant \lambda _{\max }\\r(y_{1})&\geqslant r(y_{2})\geqslant \cdots \geqslant \lambda _{\min }\end{aligned}}}$

Surge entonces la pregunta de cómo elegir los subespacios para que estas secuencias converjan a una tasa óptima.

Desde , la dirección óptima en la que buscar valores mayores de es la del gradiente , y asimismo desde la dirección óptima en la que buscar valores menores de es la del gradiente negativo . En general ${\displaystyle x_{j}}$ ${\displaystyle r}$ ${\displaystyle \nabla r(x_{j})}$ ${\displaystyle y_{j}}$ ${\displaystyle r}$ ${\displaystyle -\nabla r(y_{j})}$

${\displaystyle \nabla r(x)={\frac {2}{x^{*}x}}(Ax-r(x)x),}$

por lo tanto, las direcciones de interés son bastante fáciles de calcular en aritmética matricial, pero si uno desea mejorar en ambas y entonces hay dos nuevas direcciones a tener en cuenta: y dado que y pueden ser vectores linealmente independientes (de hecho, son casi ortogonales), en general no se puede esperar que y sean paralelos. No es necesario aumentar la dimensión de por en cada paso si se toman como subespacios de Krylov, porque entonces para todos por lo tanto en particular para ambos y . ${\displaystyle x_{j}}$ ${\displaystyle y_{j}}$ ${\displaystyle Ax_{j}}$ ${\displaystyle Ay_{j};}$ ${\displaystyle x_{j}}$ ${\displaystyle y_{j}}$ ${\displaystyle Ax_{j}}$ ${\displaystyle Ay_{j}}$ ${\displaystyle {\mathcal {L}}_{j}}$ ${\displaystyle 2}$ ${\displaystyle \{{\mathcal {L}}_{j}\}_{j=1}^{m}}$ ${\displaystyle Az\in {\mathcal {L}}_{j+1}}$ ${\displaystyle z\in {\mathcal {L}}_{j},}$ ${\displaystyle z=x_{j}}$ ${\displaystyle z=y_{j}}$

En otras palabras, podemos empezar con algún vector inicial arbitrario para construir los espacios vectoriales. ${\displaystyle x_{1}=y_{1},}$

${\displaystyle {\mathcal {L}}_{j}=\operatorname {span} (x_{1},Ax_{1},\ldots ,A^{j-1}x_{1})}$

y luego buscar tal que ${\displaystyle x_{j},y_{j}\in {\mathcal {L}}_{j}}$

${\displaystyle r(x_{j})=\max _{z\in {\mathcal {L}}_{j}}r(z)\qquad {\text{and}}\qquad r(y_{j})=\min _{z\in {\mathcal {L}}_{j}}r(z).}$

Dado que el método de potencia th pertenece a iteración , se deduce que una iteración para producir y no puede converger más lentamente que la del método de potencia, y logrará más al aproximar ambos extremos de valores propios. Para el subproblema de optimizar en algún , es conveniente tener una base ortonormal para este espacio vectorial. Por lo tanto, nuevamente nos vemos llevados al problema de calcular iterativamente dicha base para la secuencia de subespacios de Krylov. ${\displaystyle j}$ ${\displaystyle u_{j}}$ ${\displaystyle {\mathcal {L}}_{j},}$ ${\displaystyle x_{j}}$ ${\displaystyle y_{j}}$ ${\displaystyle r}$ ${\displaystyle {\mathcal {L}}_{j}}$ ${\displaystyle \{v_{1},\ldots ,v_{j}\}}$

Convergencia y otras dinámicas

Al analizar la dinámica del algoritmo, es conveniente tomar los valores propios y vectores propios de como dados, aunque no sean conocidos explícitamente por el usuario. Para fijar la notación, sean los valores propios (se sabe que todos son reales y, por lo tanto, es posible ordenarlos) y sea un conjunto ortonormal de vectores propios tales que para todo . ${\displaystyle A}$ ${\displaystyle \lambda _{1}\geqslant \lambda _{2}\geqslant \dotsb \geqslant \lambda _{n}}$ ${\displaystyle z_{1},\dotsc ,z_{n}}$ ${\displaystyle Az_{k}=\lambda _{k}z_{k}}$ ${\displaystyle k=1,\dotsc ,n}$

También es conveniente fijar una notación para los coeficientes del vector de Lanczos inicial con respecto a esta base propia; sea para todo , de modo que . Un vector de partida empobrecido de algún componente propio retrasará la convergencia al valor propio correspondiente, y aunque esto simplemente salga como un factor constante en los límites de error, el empobrecimiento sigue siendo indeseable. Una técnica común para evitar ser golpeado constantemente por él es elegir primero extrayendo los elementos aleatoriamente de acuerdo con la misma distribución normal con media y luego reescalar el vector a la norma . Antes del reescalado, esto hace que los coeficientes también sean variables estocásticas distribuidas normalmente independientes de la misma distribución normal (ya que el cambio de coordenadas es unitario), y después del reescalado el vector tendrá una distribución uniforme en la esfera unitaria en . Esto hace posible limitar la probabilidad de que, por ejemplo , . ${\displaystyle v_{1}}$ ${\displaystyle d_{k}=z_{k}^{*}v_{1}}$ ${\displaystyle k=1,\dotsc ,n}$ ${\displaystyle \textstyle v_{1}=\sum _{k=1}^{n}d_{k}z_{k}}$ ${\displaystyle v_{1}}$ ${\displaystyle v_{1}}$ ${\displaystyle 0}$ ${\displaystyle 1}$ ${\displaystyle d_{k}}$ ${\displaystyle (d_{1},\dotsc ,d_{n})}$ ${\displaystyle \mathbb {C} ^{n}}$ ${\displaystyle |d_{1}|<\varepsilon }$

El hecho de que el algoritmo de Lanczos sea independiente de las coordenadas (las operaciones solo miran los productos internos de los vectores, nunca los elementos individuales de los vectores) hace que sea fácil construir ejemplos con una estructura propia conocida para ejecutar el algoritmo: crear una matriz diagonal con los valores propios deseados en la diagonal; siempre que el vector inicial tenga suficientes elementos distintos de cero, el algoritmo generará una matriz simétrica tridiagonal general como . ${\displaystyle A}$ ${\displaystyle v_{1}}$ ${\displaystyle T}$

Teoría de convergencia de Kaniel-Paige

Después de los pasos de iteración del algoritmo de Lanczos, es una matriz simétrica real, que de manera similar a la anterior tiene valores propios. Por convergencia se entiende principalmente la convergencia de a (y la convergencia simétrica de a ) a medida que crece, y en segundo lugar la convergencia de un cierto rango de valores propios de a sus contrapartes de . La convergencia para el algoritmo de Lanczos es a menudo órdenes de magnitud más rápida que para el algoritmo de iteración de potencia. ^{[9] : 477} ${\displaystyle m}$ ${\displaystyle T}$ ${\displaystyle m\times m}$ ${\displaystyle m}$ ${\displaystyle \theta _{1}\geqslant \theta _{2}\geqslant \dots \geqslant \theta _{m}.}$ ${\displaystyle \theta _{1}}$ ${\displaystyle \lambda _{1}}$ ${\displaystyle \theta _{m}}$ ${\displaystyle \lambda _{n}}$ ${\displaystyle m}$ ${\displaystyle \theta _{1},\ldots ,\theta _{k}}$ ${\displaystyle T}$ ${\displaystyle \lambda _{1},\ldots ,\lambda _{k}}$ ${\displaystyle A}$

Los límites para provienen de la interpretación anterior de los valores propios como valores extremos del cociente de Rayleigh . Dado que es a priori el máximo de en la totalidad de mientras que es simplemente el máximo en un subespacio de Krylov de dimensión , obtenemos trivialmente . Por el contrario, cualquier punto en ese subespacio de Krylov proporciona un límite inferior para , por lo que si se puede exhibir un punto para que es pequeño, entonces esto proporciona un límite estricto en . ${\displaystyle \theta _{1}}$ ${\displaystyle r(x)}$ ${\displaystyle \lambda _{1}}$ ${\displaystyle r}$ ${\displaystyle \mathbb {C} ^{n},}$ ${\displaystyle \theta _{1}}$ ${\displaystyle m}$ ${\displaystyle \lambda _{1}\geqslant \theta _{1}}$ ${\displaystyle x}$ ${\displaystyle r(x)}$ ${\displaystyle \theta _{1}}$ ${\displaystyle \lambda _{1}-r(x)}$ ${\displaystyle \theta _{1}}$

La dimensión del subespacio de Krylov es ${\displaystyle m}$

${\displaystyle \operatorname {span} \left\{v_{1},Av_{1},A^{2}v_{1},\ldots ,A^{m-1}v_{1}\right\},}$

De modo que cualquier elemento de éste puede expresarse como para algún polinomio de grado como máximo ; los coeficientes de ese polinomio son simplemente los coeficientes en la combinación lineal de los vectores . El polinomio que queremos resultará tener coeficientes reales, pero por el momento deberíamos permitir también coeficientes complejos, y escribiremos para el polinomio obtenido mediante la conjugación compleja de todos los coeficientes de . En esta parametrización del subespacio de Krylov, tenemos ${\displaystyle p(A)v_{1}}$ ${\displaystyle p}$ ${\displaystyle m-1}$ ${\displaystyle v_{1},Av_{1},A^{2}v_{1},\ldots ,A^{m-1}v_{1}}$ ${\displaystyle p^{*}}$ ${\displaystyle p}$

${\displaystyle r(p(A)v_{1})={\frac {(p(A)v_{1})^{*}Ap(A)v_{1}}{(p(A)v_{1})^{*}p(A)v_{1}}}={\frac {v_{1}^{*}p(A)^{*}Ap(A)v_{1}}{v_{1}^{*}p(A)^{*}p(A)v_{1}}}={\frac {v_{1}^{*}p^{*}(A^{*})Ap(A)v_{1}}{v_{1}^{*}p^{*}(A^{*})p(A)v_{1}}}={\frac {v_{1}^{*}p^{*}(A)Ap(A)v_{1}}{v_{1}^{*}p^{*}(A)p(A)v_{1}}}}$

Usando ahora la expresión para como una combinación lineal de vectores propios, obtenemos ${\displaystyle v_{1}}$

${\displaystyle Av_{1}=A\sum _{k=1}^{n}d_{k}z_{k}=\sum _{k=1}^{n}d_{k}\lambda _{k}z_{k}}$

y de manera más general

${\displaystyle q(A)v_{1}=\sum _{k=1}^{n}d_{k}q(\lambda _{k})z_{k}}$

para cualquier polinomio . ${\displaystyle q}$

De este modo

${\displaystyle \lambda _{1}-r(p(A)v_{1})=\lambda _{1}-{\frac {v_{1}^{*}\sum _{k=1}^{n}d_{k}p^{*}(\lambda _{k})\lambda _{k}p(\lambda _{k})z_{k}}{v_{1}^{*}\sum _{k=1}^{n}d_{k}p^{*}(\lambda _{k})p(\lambda _{k})z_{k}}}=\lambda _{1}-{\frac {\sum _{k=1}^{n}|d_{k}|^{2}\lambda _{k}p(\lambda _{k})^{*}p(\lambda _{k})}{\sum _{k=1}^{n}|d_{k}|^{2}p(\lambda _{k})^{*}p(\lambda _{k})}}={\frac {\sum _{k=1}^{n}|d_{k}|^{2}(\lambda _{1}-\lambda _{k})\left|p(\lambda _{k})\right|^{2}}{\sum _{k=1}^{n}|d_{k}|^{2}\left|p(\lambda _{k})\right|^{2}}}.}$

Una diferencia clave entre numerador y denominador es que el término se anula en el numerador, pero no en el denominador. Por lo tanto, si se puede elegir que sea grande en pero pequeño en todos los demás valores propios, se obtendrá un límite estricto para el error . ${\displaystyle k=1}$ ${\displaystyle p}$ ${\displaystyle \lambda _{1}}$ ${\displaystyle \lambda _{1}-\theta _{1}}$

Dado que tiene muchos más valores propios que coeficientes, esto puede parecer una tarea difícil, pero una forma de cumplirla es usar polinomios de Chebyshev . Escribiendo para el grado del polinomio de Chebyshev de primera especie (el que satisface para todos ), tenemos un polinomio que permanece en el rango del intervalo conocido pero crece rápidamente fuera de él. Con algún escalamiento del argumento, podemos hacer que mapee todos los valores propios excepto en . Sea ${\displaystyle A}$ ${\displaystyle p}$ ${\displaystyle c_{k}}$ ${\displaystyle k}$ ${\displaystyle c_{k}(\cos x)=\cos(kx)}$ ${\displaystyle x}$ ${\displaystyle [-1,1]}$ ${\displaystyle [-1,1]}$ ${\displaystyle \lambda _{1}}$ ${\displaystyle [-1,1]}$

${\displaystyle p(x)=c_{m-1}\left({\frac {2x-\lambda _{2}-\lambda _{n}}{\lambda _{2}-\lambda _{n}}}\right)}$

(en caso de , utilice en su lugar el valor propio más grande estrictamente menor que ), entonces el valor máximo de para es y el valor mínimo es , por lo que ${\displaystyle \lambda _{2}=\lambda _{1}}$ ${\displaystyle \lambda _{1}}$ ${\displaystyle |p(\lambda _{k})|^{2}}$ ${\displaystyle k\geqslant 2}$ ${\displaystyle 1}$ ${\displaystyle 0}$

${\displaystyle \lambda _{1}-\theta _{1}\leqslant \lambda _{1}-r(p(A)v_{1})={\frac {\sum _{k=2}^{n}|d_{k}|^{2}(\lambda _{1}-\lambda _{k})|p(\lambda _{k})|^{2}}{\sum _{k=1}^{n}|d_{k}|^{2}|p(\lambda _{k})|^{2}}}\leqslant {\frac {\sum _{k=2}^{n}|d_{k}|^{2}(\lambda _{1}-\lambda _{k})}{|d_{1}|^{2}|p(\lambda _{1})|^{2}}}\leqslant {\frac {(\lambda _{1}-\lambda _{n})\sum _{k=2}^{n}|d_{k}|^{2}}{|p(\lambda _{1})|^{2}|d_{1}|^{2}}}.}$

Además

${\displaystyle p(\lambda _{1})=c_{m-1}\left({\frac {2\lambda _{1}-\lambda _{2}-\lambda _{n}}{\lambda _{2}-\lambda _{n}}}\right)=c_{m-1}\left(2{\frac {\lambda _{1}-\lambda _{2}}{\lambda _{2}-\lambda _{n}}}+1\right);}$

La cantidad

${\displaystyle \rho ={\frac {\lambda _{1}-\lambda _{2}}{\lambda _{2}-\lambda _{n}}}}$

(es decir, la relación entre el primer intervalo propio y el diámetro del resto del espectro ) es, por lo tanto, de importancia clave para la tasa de convergencia aquí. También se escribe

${\displaystyle R=e^{\operatorname {arcosh} (1+2\rho )}=1+2\rho +2{\sqrt {\rho ^{2}+\rho }},}$

Podemos concluir que

${\displaystyle {\begin{aligned}\lambda _{1}-\theta _{1}&\leqslant {\frac {(\lambda _{1}-\lambda _{n})\left(1-|d_{1}|^{2}\right)}{c_{m-1}(2\rho +1)^{2}|d_{1}|^{2}}}\\[6pt]&={\frac {1-|d_{1}|^{2}}{|d_{1}|^{2}}}(\lambda _{1}-\lambda _{n}){\frac {1}{\cosh ^{2}((m-1)\operatorname {arcosh} (1+2\rho ))}}\\[6pt]&={\frac {1-|d_{1}|^{2}}{|d_{1}|^{2}}}(\lambda _{1}-\lambda _{n}){\frac {4}{\left(R^{m-1}+R^{-(m-1)}\right)^{2}}}\\[6pt]&\leqslant 4{\frac {1-|d_{1}|^{2}}{|d_{1}|^{2}}}(\lambda _{1}-\lambda _{n})R^{-2(m-1)}\end{aligned}}}$

De este modo, la tasa de convergencia está controlada principalmente por , ya que este límite se reduce en un factor por cada iteración adicional. ${\displaystyle R}$ ${\displaystyle R^{-2}}$

A modo de comparación, se puede considerar cómo la tasa de convergencia del método de potencia depende de , pero dado que el método de potencia es principalmente sensible al cociente entre los valores absolutos de los valores propios, necesitamos que la brecha propia entre y sea la dominante. Bajo esa restricción, el caso que más favorece al método de potencia es que , así que considere que. Más adelante en el método de potencia, el vector de iteración: ${\displaystyle \rho }$ ${\displaystyle |\lambda _{n}|\leqslant |\lambda _{2}|}$ ${\displaystyle \lambda _{1}}$ ${\displaystyle \lambda _{2}}$ ${\displaystyle \lambda _{n}=-\lambda _{2}}$

${\displaystyle u=(1-t^{2})^{1/2}z_{1}+tz_{2}\approx z_{1}+tz_{2},}$^{[nota 1]}

donde cada nueva iteración multiplica efectivamente la amplitud por ${\displaystyle z_{2}}$ ${\displaystyle t}$

${\displaystyle {\frac {\lambda _{2}}{\lambda _{1}}}={\frac {\lambda _{2}}{\lambda _{2}+(\lambda _{1}-\lambda _{2})}}={\frac {1}{1+{\frac {\lambda _{1}-\lambda _{2}}{\lambda _{2}}}}}={\frac {1}{1+2\rho }}.}$

La estimación del valor propio más grande es entonces

${\displaystyle u^{*}Au=(1-t^{2})\lambda _{1}+t^{2}\lambda _{2},}$

Por lo tanto, el límite anterior para la tasa de convergencia del algoritmo de Lanczos debe compararse con

${\displaystyle \lambda _{1}-u^{*}Au=(\lambda _{1}-\lambda _{2})t^{2},}$

que se reduce en un factor de para cada iteración. La diferencia, por tanto, se reduce a la que existe entre y . En la región, este último es más parecido a , y funciona como lo haría el método de potencia con un eigengap dos veces más grande; una mejora notable. Sin embargo, el caso más desafiante es el de en el que hay una mejora aún mayor en el eigengap; la región es donde el algoritmo de Lanczos en términos de convergencia hace la mejora más pequeña en el método de potencia. ${\displaystyle (1+2\rho )^{-2}}$ ${\displaystyle 1+2\rho }$ ${\displaystyle R=1+2\rho +2{\sqrt {\rho ^{2}+\rho }}}$ ${\displaystyle \rho \gg 1}$ ${\displaystyle 1+4\rho }$ ${\displaystyle \rho \ll 1,}$ ${\displaystyle R\approx 1+2{\sqrt {\rho }}}$ ${\displaystyle \rho \gg 1}$

Estabilidad numérica

La estabilidad se refiere a cuánto se verá afectado el algoritmo (es decir, si producirá un resultado aproximado cercano al original) si se introducen y acumulan pequeños errores numéricos. La estabilidad numérica es el criterio central para juzgar la utilidad de implementar un algoritmo en una computadora con redondeo.

En el caso del algoritmo de Lanczos, se puede demostrar que con aritmética exacta , el conjunto de vectores construye una base ortonormal , y los valores propios/vectores resueltos son buenas aproximaciones a los de la matriz original. Sin embargo, en la práctica (ya que los cálculos se realizan en aritmética de punto flotante donde la inexactitud es inevitable), la ortogonalidad se pierde rápidamente y en algunos casos el nuevo vector podría incluso ser linealmente dependiente del conjunto que ya está construido. Como resultado, algunos de los valores propios de la matriz tridiagonal resultante pueden no ser aproximaciones a la matriz original. Por lo tanto, el algoritmo de Lanczos no es muy estable. ${\displaystyle v_{1},v_{2},\cdots ,v_{m+1}}$

Los usuarios de este algoritmo deben poder encontrar y eliminar esos valores propios "espurios". Las implementaciones prácticas del algoritmo de Lanczos van en tres direcciones para combatir este problema de estabilidad: ^{[6] [7]}

1. Evitar la pérdida de ortogonalidad,
2. Recuperar la ortogonalidad después de generar la base.
3. Una vez identificados todos los valores propios buenos y "espurios", elimine los espurios.

Variaciones

Existen variaciones del algoritmo de Lanczos en las que los vectores implicados son matrices altas y estrechas en lugar de vectores y las constantes de normalización son matrices cuadradas pequeñas. Estos algoritmos se denominan algoritmos de Lanczos de "bloque" y pueden ser mucho más rápidos en computadoras con una gran cantidad de registros y tiempos de recuperación de memoria prolongados.

Muchas implementaciones del algoritmo Lanczos se reinician después de una cierta cantidad de iteraciones. Una de las variaciones reiniciadas más influyentes es el método Lanczos reiniciado implícitamente, ^[10] que se implementa en ARPACK . ^[11] Esto ha dado lugar a una serie de otras variaciones reiniciadas, como la bidiagonalización Lanczos reiniciada. ^[12] Otra variación reiniciada exitosa es el método Lanczos Thick-Restart, ^[13] que se ha implementado en un paquete de software llamado TRLan. ^[14]

Espacio nulo sobre un cuerpo finito

En 1995, Peter Montgomery publicó un algoritmo, basado en el algoritmo de Lanczos, para encontrar elementos del espacio nulo de una matriz dispersa grande sobre GF(2) ; dado que el conjunto de personas interesadas en matrices dispersas grandes sobre campos finitos y el conjunto de personas interesadas en problemas de valores propios grandes apenas se superponen, a menudo esto también se denomina algoritmo de bloque de Lanczos sin causar una confusión irrazonable. ^{[ cita requerida ]}

Aplicaciones

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

Los algoritmos de Lanczos también se utilizan en física de materia condensada como método para resolver hamiltonianos de sistemas electrónicos fuertemente correlacionados , ^[15] así como en códigos de modelos de capas en física nuclear . ^[16]

Implementaciones

La biblioteca NAG contiene varias rutinas ^[17] para la solución de sistemas lineales de gran escala y problemas propios que utilizan el algoritmo de Lanczos.

MATLAB y GNU Octave vienen con ARPACK incorporado. Tanto las matrices almacenadas como las implícitas se pueden analizar mediante la función eigs() (Matlab/Octave).

De manera similar, en Python , el paquete SciPy tiene scipy.sparse.linalg.eigsh, que también es un contenedor para las funciones SSEUPD y DSEUPD de ARPACK que utilizan el método Lanczos reiniciado implícitamente.

Una implementación de Matlab del algoritmo Lanczos (nótese los problemas de precisión) está disponible como parte del paquete de Matlab de propagación de creencias gaussianas. La biblioteca de filtrado colaborativo GraphLab ^[18] incorpora una implementación paralela a gran escala del algoritmo Lanczos (en C++) para núcleos múltiples.

La biblioteca PRIMME también implementa un algoritmo similar a Lanczos.

Notas

1. No es necesario que ambos coeficientes sean reales, pero la fase tiene poca importancia. Tampoco es necesario que los componentes de otros vectores propios hayan desaparecido por completo, pero se reducen al menos tan rápido como el de , por lo que describe el peor caso. ${\displaystyle z_{2}}$ ${\displaystyle u\approx z_{1}+tz_{2}}$

Referencias

1. Lanczos, C. (1950). "Un método de iteración para la solución del problema de valores propios de operadores diferenciales e integrales lineales" (PDF) . Revista de investigación de la Oficina Nacional de Normas . 45 (4): 255–282. doi : 10.6028/jres.045.026 .
2. Ojalvo, IU; Newman, M. (1970). "Modos de vibración de grandes estructuras mediante un método automático de reducción matricial". AIAA Journal . 8 (7): 1234–1239. Bibcode :1970AIAAJ...8.1234N. doi :10.2514/3.5878.
3. Paige, CC (1971). El cálculo de valores propios y vectores propios de matrices dispersas muy grandes (tesis doctoral). Universidad de Londres. OCLC  654214109.
4. Paige, CC (1972). "Variantes computacionales del método de Lanczos para el problema propio". J. Inst. Maths Applics . 10 (3): 373–381. doi :10.1093/imamat/10.3.373.
5. Ojalvo, IU (1988). "Orígenes y ventajas de los vectores de Lanczos para grandes sistemas dinámicos". Proc. 6th Modal Analysis Conference (IMAC), Kissimmee, FL . págs. 489–494.
6. Cullum; Willoughby (1985). Algoritmos de Lanczos para cálculos de valores propios simétricos grandes . Vol. 1. ISBN 0-8176-3058-9.
7. Yousef Saad (22 de junio de 1992). Métodos numéricos para problemas de valores propios grandes. ISBN 0-470-21820-7.
8. Coakley, Ed S.; Rokhlin, Vladimir (2013). "Un algoritmo rápido de divide y vencerás para calcular los espectros de matrices tridiagonales simétricas reales". Análisis armónico computacional y aplicado . 34 (3): 379–414. doi :10.1016/j.acha.2012.06.003.
9. Golub, Gene H.; Van Loan, Charles F. (1996). Cálculos matriciales (3.ª ed.). Baltimore: Johns Hopkins Univ. Press. ISBN 0-8018-5413-X.
10. D. Calvetti ; L. Reichel; DC Sorensen (1994). "Un método de Lanczos reiniciado implícitamente para problemas de valores propios simétricos grandes". Transacciones electrónicas sobre análisis numérico . 2 : 1–21.
11. RB Lehoucq; DC Sorensen; C. Yang (1998). Guía del usuario de ARPACK: Solución de problemas de valores propios a gran escala con métodos de Arnoldi reiniciados implícitamente . SIAM. doi :10.1137/1.9780898719628. ISBN . 978-0-89871-407-4.
12. E. Kokiopoulou; C. Bekas; E. Gallopoulos (2004). "Cálculo de los tripletes singulares más pequeños con bidiagonalización de Lanczos reiniciada implícitamente" (PDF) . Appl. Numer. Math . 49 : 39–61. doi :10.1016/j.apnum.2003.11.011.
13. Kesheng Wu; Horst Simon (2000). "Método de Lanczos de reinicio grueso para problemas de valores propios simétricos grandes". Revista SIAM sobre análisis de matrices y aplicaciones . 22 (2). SIAM: 602–616. doi :10.1137/S0895479898334605.
14. Kesheng Wu; Horst Simon (2001). «Paquete de software TRLan». Archivado desde el original el 1 de julio de 2007. Consultado el 30 de junio de 2007 .
15. Chen, HY; Atkinson, WA; Wortis, R. (julio de 2011). "Anomalía de sesgo cero inducida por desorden en el modelo de Anderson-Hubbard: cálculos numéricos y analíticos". Physical Review B . 84 (4): 045113. arXiv : 1012.1031 . Bibcode :2011PhRvB..84d5113C. doi :10.1103/PhysRevB.84.045113. S2CID  118722138.
16. Shimizu, Noritaka (21 de octubre de 2013). "Código de modelo de capas nucleares para computación paralela masiva, "KSHELL"". arXiv : 1310.5431 [nucl-th].
17. The Numerical Algorithms Group. «Índice de palabras clave: Lanczos». Manual de la biblioteca NAG, Mark 23. Consultado el 9 de febrero de 2012 .
18. GraphLab Archivado el 14 de marzo de 2011 en Wayback Machine.

Lectura adicional

• Golub, Gene H. ; Van Loan, Charles F. (1996). "Métodos Lanczos". Cálculos matriciales . Baltimore: Johns Hopkins University Press. págs. 470–507. ISBN 0-8018-5414-8.
• Ng, Andrew Y .; Zheng, Alice X.; Jordan, Michael I. (2001). "Análisis de enlaces, vectores propios y estabilidad" (PDF) . IJCAI'01 Actas de la 17.ª Conferencia conjunta internacional sobre inteligencia artificial . 2 : 903–910.
• Erik Koch (2019). «Diagonalización Exacta y Método Lanczos» (PDF) . En E. Pavarini; E. Koch; S. Zhang (eds.). "Métodos de muchos cuerpos para materiales reales" . Jülich. ISBN 978-3-95806-400-3.