Algoritmo de Lanczos

El algoritmo de Lanczos se suele mencionar para hallar los autovalores y autovectores de una matriz. Sin embargo, mientras que una diagonalización convencional de una matriz permite identificar los autovalores y autovectores a simple vista, esto no ocurre con la tridiagonalización que realiza el algoritmo de Lanczos; se requieren pasos adicionales complejos para calcular incluso un solo autovalor o autovector. No obstante, la aplicación del algoritmo de Lanczos suele representar un avance significativo en el cálculo de la descomposición en autovalores.

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

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

De este modo, el algoritmo de Lanczos transforma el problema de descomposición en valores propios para${\displaystyle A}$en el problema de descomposición de valores propios para${\displaystyle T}$.

1. Para matrices tridiagonales, existen varios algoritmos especializados, a menudo con una complejidad computacional mejor que la de los algoritmos de propósito general. Por ejemplo, si${\displaystyle T}$es un${\displaystyle m\times m}$Matriz simétrica tridiagonal entonces:
   • La recursión continua permite calcular el polinomio característico en${\displaystyle O(m^{2})}$operaciones y evaluarlo en un punto en${\displaystyle O(m)}$operaciones.
   • El algoritmo de autovalores de divide y vencerás se puede utilizar para calcular la descomposición en autovalores completa de${\displaystyle T}$en${\displaystyle O(m^{2})}$operaciones.
   • El método multipolar rápido ^{[ 8 ]} puede calcular todos los valores propios en tan solo${\displaystyle O(m\log m)}$operaciones.
2. Se sabe que algunos algoritmos generales de descomposición en valores 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${\displaystyle O(m^{2})}$igual que para el algoritmo de divide y vencerás (aunque el factor constante puede ser diferente); puesto que los autovectores juntos tienen${\displaystyle m^{2}}$elementos, esto es asintóticamente óptimo .
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 obtener beneficios de rendimiento de bajo nivel al aplicarse a la matriz tridiagonal.${\displaystyle T}$en lugar de la matriz original${\displaystyle A}$. Desde${\displaystyle T}$es muy disperso 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é . Asimismo,${\displaystyle T}$es una matriz real con todos los vectores propios y valores propios reales, mientras que${\displaystyle A}$En general, puede tener elementos y vectores propios complejos, por lo que la aritmética real es suficiente para encontrar los vectores propios y los valores propios de${\displaystyle T}$.
4. Si${\displaystyle n}$es muy grande, luego reducir${\displaystyle m}$de modo que${\displaystyle T}$es de un tamaño manejable aún permitirá encontrar los valores propios y vectores propios más extremos de${\displaystyle A}$; en el${\displaystyle m\ll n}$En esta región, el algoritmo de Lanczos puede considerarse un esquema de compresión con pérdidas para matrices hermíticas, que hace hincapié en la preservación de los valores propios extremos.

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

Aunque el problema de valores propios suele ser la motivación para aplicar el algoritmo de Lanczos, la operación que realiza principalmente es la tridiagonalización de una matriz, para la cual se han preferido las transformaciones de Householder numéricamente estables desde la década de 1950. Durante la década de 1960, el algoritmo de Lanczos cayó en desuso. El interés en él se reavivó con 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 la alternativa que se prueba solo si Householder no resulta satisfactorio. ^{[ 9 ]}

Entre los aspectos en los que difieren los dos algoritmos se incluyen:

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

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

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

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

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.${\displaystyle A}$pero solo presta atención al último resultado; las implementaciones suelen usar la misma variable para todos los vectores.${\displaystyle u_{j}}$En este caso, cada nueva iteración sobrescribe los resultados de la anterior. Quizás sea preferible conservar todos los resultados intermedios y organizar los datos.

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

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

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

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

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

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

${\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 el${\displaystyle u_{j}}$vectores para calcular estos${\displaystyle v_{j}}$, porque${\displaystyle u_{j}-v_{j}\in \operatorname {span} (v_{1},\dotsc ,v_{j-1})}$y por lo tanto la diferencia entre${\displaystyle u_{j+1}'=Au_{j}}$y${\displaystyle w_{j+1}'=Av_{j}}$está en${\displaystyle \operatorname {span} (v_{1},\dotsc ,v_{j})}$, 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

1. Elige un vector aleatorio${\displaystyle v_{1}}$de norma euclidiana${\displaystyle 1}$.
2. Para${\displaystyle j=1,\dotsc ,m-1}$hacer:
   1. Dejar${\displaystyle w_{j+1}'=Av_{j}}$.
   2. A pesar de${\displaystyle k=1,\dotsc ,j}$dejar${\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${\displaystyle h_{j+1,j}\neq 0}$entonces deja${\displaystyle v_{j+1}=w_{j+1}/h_{j+1,j}}$,

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

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

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

la definición${\displaystyle h_{j+1,j}=\|w_{j+1}\|}$Puede parecer un poco extraño, pero encaja con el patrón general.${\displaystyle h_{k,j}=v_{k}^{*}w_{j+1}'}$desde

${\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 los vectores de iteración de potencia${\displaystyle u_{j}}$que fueron eliminados de esta recursión satisfacen${\displaystyle u_{j}\in \operatorname {span} (v_{1},\ldots ,v_{j}),}$los vectores${\displaystyle \{v_{j}\}_{j=1}^{m}}$y coeficientes${\displaystyle h_{k,j}}$contienen suficiente información de${\displaystyle A}$que todo${\displaystyle u_{1},\ldots ,u_{m}}$Se puede calcular, por lo que no se perdió nada al cambiar los vectores. (De hecho, resulta que los datos recopilados aquí proporcionan aproximaciones significativamente mejores del mayor valor propio que las que se obtienen con un número igual de iteraciones en el método de potencia, aunque esto no sea necesariamente obvio en este punto).

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

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

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

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

${\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. Para${\displaystyle k=j}$uno consigue

${\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 forma más abstracta, si${\displaystyle V}$es la matriz con columnas${\displaystyle v_{1},\ldots ,v_{m}}$luego los números${\displaystyle h_{k,j}}$pueden identificarse como elementos de la matriz${\displaystyle H=V^{*}AV}$, y${\displaystyle h_{k,j}=0}$para${\displaystyle k>j+1;}$la matriz${\displaystyle H}$es el Hessenberg superior . Desde

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

la matriz${\displaystyle H}$es hermitiano. Esto implica que${\displaystyle H}$También es Hessenberg inferior, por lo que de hecho debe ser tridiagnática. Al ser hermitiana, su diagonal principal es real, y dado que su primera subdiagonal es real por construcción, lo mismo ocurre con su primera superdiagonal. Por lo tanto,${\displaystyle H}$es una matriz real y simétrica: la matriz${\displaystyle T}$de la especificación del algoritmo de Lanczos.

Una forma de caracterizar los autovectores de una matriz hermitiana${\displaystyle A}$es como puntos estacionarios del cociente de Rayleigh

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

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

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

${\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 de manera que estas secuencias converjan a una velocidad óptima.

De${\displaystyle x_{j}}$, la dirección óptima en la que buscar valores mayores de${\displaystyle r}$es la del gradiente${\displaystyle \nabla r(x_{j})}$y asimismo de${\displaystyle y_{j}}$la dirección óptima en la que buscar valores más pequeños de${\displaystyle r}$es la del gradiente negativo${\displaystyle -\nabla r(y_{j})}$. En general

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

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

En otras palabras, podemos comenzar con un vector inicial arbitrario.${\displaystyle x_{1}=y_{1},}$construir los espacios vectoriales

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

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

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

Desde el${\displaystyle j}$método de potencia iterar${\displaystyle u_{j}}$pertenece a${\displaystyle {\mathcal {L}}_{j},}$De ello se deduce que una iteración para producir el${\displaystyle x_{j}}$y${\displaystyle y_{j}}$no puede converger más lentamente que el método de potencia y logrará más aproximando ambos extremos de los valores propios. Para el subproblema de optimización${\displaystyle r}$en algunos${\displaystyle {\mathcal {L}}_{j}}$Es conveniente tener una base ortonormal.${\displaystyle \{v_{1},\ldots ,v_{j}\}}$para este espacio vectorial . Así, volvemos al problema de calcular iterativamente dicha base para la secuencia de subespacios de Krylov.

Después${\displaystyle m}$pasos de iteración del algoritmo de Lanczos,${\displaystyle T}$es un${\displaystyle m\times m}$matriz simétrica real, que de manera similar a la anterior tiene${\displaystyle m}$valores propios${\displaystyle \theta _{1}\geqslant \theta _{2}\geqslant \dots \geqslant \theta _{m}.}$La convergencia se entiende principalmente como convergencia de${\displaystyle \theta _{1}}$a${\displaystyle \lambda _{1}}$(y la convergencia simétrica de${\displaystyle \theta _{m}}$a${\displaystyle \lambda _{n}}$) como${\displaystyle m}$crece y, secundariamente, la convergencia de algún rango${\displaystyle \theta _{1},\ldots ,\theta _{k}}$de valores propios de${\displaystyle T}$a sus homólogos${\displaystyle \lambda _{1},\ldots ,\lambda _{k}}$de${\displaystyle A}$La convergencia del algoritmo de Lanczos suele ser órdenes de magnitud más rápida que la del algoritmo de iteración de potencias. ^{[ 9 ] : 477}

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

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

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

por lo que cualquier elemento del mismo puede expresarse como${\displaystyle p(A)v_{1}}$para algún polinomio${\displaystyle p}$de grado como máximo${\displaystyle m-1}$; los coeficientes de ese polinomio son simplemente los coeficientes de la combinación lineal de los vectores${\displaystyle v_{1},Av_{1},A^{2}v_{1},\ldots ,A^{m-1}v_{1}}$El polinomio que buscamos tendrá coeficientes reales, pero por el momento debemos considerar también coeficientes complejos, y escribiremos: .${\displaystyle p^{*}}$para el polinomio obtenido al conjugar de forma compleja todos los coeficientes de${\displaystyle p}$. En esta parametrización del subespacio de Krylov, tenemos

${\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${\displaystyle v_{1}}$como una combinación lineal de autovectores, obtenemos

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

y de forma 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 aquí es que${\displaystyle k=1}$El término desaparece en el numerador, pero no en el denominador. Por lo tanto, si se puede elegir${\displaystyle p}$ser grande en${\displaystyle \lambda _{1}}$pero pequeño en todos los demás autovalores, se obtendrá una cota ajustada para el error.${\displaystyle \lambda _{1}-\theta _{1}}$.

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

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

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

${\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 eigengap y el diámetro del resto del espectro ) es, por lo tanto, de vital importancia para la tasa de convergencia en este caso. También escribir

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

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

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

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

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

${\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 mayor valor propio es entonces

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

por lo que 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 por un factor de${\displaystyle (1+2\rho )^{-2}}$para cada iteración. La diferencia, por lo tanto, se reduce a que entre${\displaystyle 1+2\rho }$y${\displaystyle R=1+2\rho +2{\sqrt {\rho ^{2}+\rho }}}$. En el${\displaystyle \rho \gg 1}$región, esta última es más parecida${\displaystyle 1+4\rho }$y se comporta como lo haría el método de potencia con una brecha propia dos veces mayor; una mejora notable. Sin embargo, el caso más desafiante es el de${\displaystyle \rho \ll 1,}$en el cual${\displaystyle R\approx 1+2{\sqrt {\rho }}}$es una mejora aún mayor en el eigengap; el${\displaystyle \rho \gg 1}$La región es donde el algoritmo de Lanczos, en términos de convergencia, realiza la menor mejora con respecto al método de potencia.

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 cuerpos finitos y el conjunto de personas interesadas en problemas de valores propios grandes apenas se superponen, a menudo también se le llama algoritmo de Lanczos por bloques sin causar una confusión irrazonable.