Algoritmo de Lanczos

format_list_bulleted Contenido keyboard_arrow_down

ImprimirCitar

Cálculo Numérico eigenvalue

El algoritmo de lanczos es un método iterativo ideado por Cornelius Lanczos que es una adaptación de métodos de poder para encontrar el ${displaystyle m}$ eigenvalues y eigenvectores de un ${displaystyle ntimes n}$ Matriz hermitiana, donde ${displaystyle m}$ es a menudo pero no necesariamente mucho más pequeño que ${displaystyle n}$ . Aunque en principio era computacionalmente eficiente, el método tal como se formuló inicialmente no era útil, debido a su inestabilidad numérica.

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 de gran tamaño sometidas a cargas dinámicas. 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) con cualquier grado de precisión, que cuando no se realizaba, producía una serie de vectores que estaban altamente contaminados por aquellos asociados con las frecuencias naturales más bajas.

En su obra original, estos autores también sugirieron cómo seleccionar un vector inicial (es decir, utilizar un generador de números aleatorios para seleccionar cada elemento del vector inicial) y sugirieron un método empíricamente determinado para determinar ${displaystyle m}$ , el número reducido de vectores (es decir, debe ser seleccionado para ser aproximadamente 1,5 veces el número de eigenvalues precisos deseados). Poco después su trabajo fue seguido por Paige, quien también proporcionó un análisis de errores. En 1988, Ojalvo produjo una historia más detallada de este algoritmo y una prueba de error de eigenvalue eficiente.

El algoritmo

Input a Hermitian matriz

{displaystyle A}

de tamaño

{displaystyle ntimes n}

, y opcionalmente varias iteraciones

{displaystyle m}

(como predeterminado, dejar

{displaystyle m=n}

Estrictamente hablando, el algoritmo no necesita acceso a la matriz explícita, pero sólo una función ${displaystyle vmapsto Av}$ que computa el producto de la matriz por un vector arbitrario. Esta función se llama a la mayoría ${displaystyle m}$ veces.

Producto an

{displaystyle ntimes m}

matriz

{displaystyle V}

con columnas ortonormales y una matriz simétrica real tridiagonal

{displaystyle T=V^{*}AV}

de tamaño

{displaystyle mtimes m}

. Si

{displaystyle m=n}

Entonces

{displaystyle V}

es unitario, y

{displaystyle A=VTV^{*}}

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

Vamos. ${displaystyle v_{1}in mathbb {C} ^{n}}$ ser un vector arbitrario con la norma Euclidea ${displaystyle 1}$ .
Paso de iteración inicial Abreviada:
1. Vamos. ${displaystyle w_{1}'=Av_{1}}$ .
2. Vamos. ${displaystyle alpha _{1}=w_{1}'^{*}v_{1}}$ .
3. Vamos. ${displaystyle w_{1}=w_{1}'-alpha _{1}v_{1}}$ .
Para j=2,... ... ,m{displaystyle j=2,dotsm} ${displaystyle j=2,dotsm}$ hacer:
1. Vamos. ${displaystyle beta _{j}=|w_{j-1}|}$ (también la norma Euclideana).
2. Si ${displaystyle beta _{j}neq 0}$ , entonces deja ${displaystyle v_{j}=w_{j-1}/beta _{j}}$ ,
  más elegir como ${displaystyle v_{j}}$ un vector arbitrario con la norma Euclidea ${displaystyle 1}$ que es ortogonal a todos ${displaystyle v_{1},dotsv_{j-1}}$ .
3. Vamos. ${displaystyle w_{j}'=Av_{j}}$ .
4. Vamos. ${displaystyle alpha _{j}=w_{j}'^{*}v_{j}}$ .
5. Vamos. ${displaystyle w_{j}=w_{j}'-alpha _{j}v_{j}-beta _{j}v_{j-1}}$ .
Vamos. ${displaystyle V}$ ser la matriz con columnas ${displaystyle v_{1},dotsv_{m}}$ . Vamos. ${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

{displaystyle Av_{j}=w_{j}'=beta _{j+1}v_{j+1}+alpha _{j}v_{j}+beta _{j}v_{j-1}}

para

<math alttext="{displaystyle 2<j2c)jc)m{displaystyle 2 segÃ3n] <img alt="{displaystyle 2<j

Hay en principio cuatro maneras de escribir el procedimiento de iteración. Paige y otros trabajos muestran que el orden de operaciones anterior es el más numérico estable. En la práctica el vector inicial ${displaystyle v_{1}}$ puede ser tomado como otro argumento del procedimiento, con ${displaystyle beta _{j}=0}$ e indicadores de imprecisión numérica que se incluyen como condiciones adicionales de terminación de bucle.

Sin contar la multiplicación de la matriz-vector, cada iteración lo hace ${displaystyle O(n)}$ operaciones aritméticas. La multiplicación de la matriz-vector se puede hacer en ${displaystyle O(dn)}$ operaciones aritméticas donde ${displaystyle d}$ es el número promedio de elementos no cero seguidos. La complejidad total es así ${displaystyle O(dmn)}$ o ${displaystyle O(dn^{2})}$ si ${displaystyle m=n}$ ; el algoritmo Lanczos puede ser muy rápido para las matrices escasas. Los esquemas para mejorar la estabilidad numérica suelen juzgarse contra este alto rendimiento.

Los vectores ${displaystyle v_{j}}$ se llaman vectores Lanczos. El vector ${displaystyle w_{j}'}$ no se utiliza después ${displaystyle w_{j}}$ es computado, y el vector ${displaystyle w_{j}}$ no se utiliza después ${displaystyle v_{j+1}}$ está computado. Por lo tanto, uno puede utilizar el mismo almacenamiento para los tres. Del mismo modo, si sólo la matriz tridiagonal ${displaystyle T}$ se busca, entonces la iteración cruda no necesita ${displaystyle v_{j-1}}$ después de haber computado ${displaystyle w_{j}}$ , aunque algunos planes para mejorar la estabilidad numérica lo necesitarían más adelante. A veces los vectores Lanczos posteriores son recomputados de ${displaystyle v_{1}}$ cuando sea necesario.

Aplicación al problema propio

El algoritmo de Lanczos se plantea con más frecuencia en el contexto de encontrar los valores eigenvectores y eigenvectores de una matriz, pero mientras que una diagonalización ordinaria de una matriz haría que los eigenvectores y los eigenvalues aparentes de la inspección, lo mismo no es cierto para la tridiagonalización realizada por el algoritmo de Lanczos; se necesitan pasos adicionales notriviales para computar un solo eige Sin embargo, aplicar el algoritmo Lanczos es a menudo un avance significativo en la computación de la eigendecomposición. Si ${displaystyle lambda }$ es un eigenvalue de ${displaystyle A}$ Y si ${displaystyle Tx=lambda x}$ () ${displaystyle x}$ es un eigenvector de ${displaystyle T}$ entonces ${displaystyle y=Vx}$ es el eigenvector correspondiente ${displaystyle A}$ (since ${displaystyle Ay=AVx=VTV^{*}Vx=VTIx=VTx=V(lambda x)=lambda Vx=lambda y}$ ). Así el algoritmo Lanczos transforma el problema de la eigendecomposición para ${displaystyle A}$ en el problema de la eigendecomposición ${displaystyle T}$ .

Para las matrices tridiagonales existen varios algoritmos especializados, a menudo con mejor complejidad computacional que algoritmos de uso general. Por ejemplo, si T{displaystyle T} ${displaystyle T}$ es un m× × m{displaystyle mtimes m} ${displaystyle mtimes m}$ matriz simétrica tridiagonal entonces:
- La recurrencia continua permite calcular el polinomio característico en ${displaystyle O(m^{2})}$ operaciones y evaluación en un punto ${displaystyle O(m)}$ operaciones.
- El algoritmo eigenvalue divide y conquista puede utilizarse para calcular toda la eigendecomposición de ${displaystyle T}$ dentro ${displaystyle O(m^{2})}$ operaciones.
- El Multipollo rápido El método puede calcular todos los eigenvalues en sólo ${displaystyle O(mlog m)}$ operaciones.
Algunos algoritmos generales de eigendecomposición, en particular el algoritmo QR, se sabe que convergen más rápido para matrices tridiagonales que para matrices generales. La complejidad asintotica de QR tridiagonal es ${displaystyle O(m^{2})}$ así como para el algoritmo de división y conquista (aunque el factor constante puede ser diferente); ya que los eigenvectores juntos tienen ${displaystyle m^{2}}$ elementos, esto es asintóticamente óptimo.
Incluso algoritmos cuyas tasas de convergencia no son afectadas por transformaciones unitarias, como el método de potencia y la iteración inversa, pueden disfrutar de beneficios de bajo nivel de rendimiento de ser aplicados a la matriz tridiagonal ${displaystyle T}$ en lugar de la matriz original ${displaystyle A}$ . Desde ${displaystyle T}$ es muy escaso con todos los elementos no cero en posiciones altamente predecibles, permite un almacenamiento compacto con excelente rendimiento frente al caché. Del mismo modo, ${displaystyle T}$ es una matriz real con todos los eigenvectores y eigenvalues reales, mientras que ${displaystyle A}$ en general puede tener elementos complejos y eigenvectores, por lo que la aritmética real es suficiente para encontrar los eigenvectores y eigenvalues de ${displaystyle T}$ .
Si ${displaystyle n}$ es muy grande, luego reducir ${displaystyle m}$ así ${displaystyle T}$ es de un tamaño manejable todavía permitirá encontrar los eigenvalues más extremos y eigenvectores de ${displaystyle A}$ ; en ${displaystyle mll n}$ región, el algoritmo Lanczos se puede ver como un esquema de compresión perdido para matrices hermitianas, que enfatiza la preservación de los eigenvalues extremos.

La combinación de buen rendimiento para matrices dispersas y la capacidad de calcular varios (sin calcular todos) los valores propios 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 realiza principalmente el algoritmo es la tridiagonalización de una matriz, para la cual las transformaciones de Householder numéricamente estables se han favorecido desde la década de 1950. Durante la década de 1960, se hizo caso omiso del algoritmo de Lanczos. El interés en él se vio rejuvenecido por la teoría de la 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 sólo si Householder no es satisfactorio.

Los aspectos en los que difieren los dos algoritmos incluyen:

Lanczos se aprovecha de ${displaystyle A}$ ser una matriz escasa, mientras que el dueño de la casa no, y generará relleno.
Lanczos trabaja con la matriz original ${displaystyle A}$ (y no tiene ningún problema con que se le conozca sólo implícitamente), mientras que el titular de la Casa cruda quiere modificar la matriz durante la computación (aunque eso puede evitarse).
Cada iteración del algoritmo Lanczos produce otra columna de la matriz de transformación final ${displaystyle V}$ , mientras que una iteración de los propietarios de la casa produce otro factor en una factorización unitaria ${displaystyle Q_{1}Q_{2}dots Q_{n}}$ de ${displaystyle V}$ . Cada factor está determinado por un solo vector, por lo que los requisitos de almacenamiento son los mismos para ambos algoritmos, y ${displaystyle V=Q_{1}Q_{2}dots Q_{n}}$ puede ser calculado en ${displaystyle O(n^{3})}$ tiempo.
El dueño de la casa es numéricamente estable, mientras que el Lanczos crudo no es.
Lanczos es muy paralelo, con sólo ${displaystyle O(n)}$ puntos de sincronización (las computaciones de ${displaystyle alpha _{j}}$ y ${displaystyle beta _{j}}$ ). El dueño de la casa es menos paralelo, teniendo una secuencia de ${displaystyle O(n^{2})}$ cantidades de escalar calculadas que cada uno depende de la cantidad anterior en la secuencia.

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 eigenvalo de mayor magnitud y un eigenvector correspondiente de una matriz ${displaystyle A}$ es difícil

Elija un vector aleatorio ${displaystyle u_{1}neq 0}$ .
Para j⩾ ⩾ 1{displaystyle jgeqslant 1} ${displaystyle jgeqslant 1}$ (hasta la dirección uj{displaystyle u_{j} ${displaystyle u_{j}}$ ha convergedo) hacer:
1. Vamos. ${displaystyle u_{j+1}'=Au_{j}.}$
2. Vamos. ${displaystyle u_{j+1}=u_{j+1}'/|u_{j+1}'|.}$

En grande ${displaystyle j}$ límite, ${displaystyle u_{j}}$ se acerca al eigenvector normado correspondiente al eigenvalu de mayor magnitud.

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

Un pedazo de información que trivialmente está disponible de los vectores ${displaystyle u_{j}}$ es una cadena de subespacios Krylov. Una manera de afirmar que sin introducir conjuntos en el algoritmo es afirmar que compute

a subconjunto

{displaystyle {v_{j}}_{j=1}^{m}}

de una base

{displaystyle mathbb {C} ^{n}}

tales que

{displaystyle Axin operatorname {span} (v_{1},dotscv_{j+1})}

para todos

{displaystyle xin operatorname {span} (v_{1},dotscv_{j})}

y todos

<math alttext="{displaystyle 1leqslant j1⩽ ⩽ jc)m;{displaystyle 1leqslant jierem;} <img alt="{displaystyle 1leqslant j

esto es trivialmente satisfecho ${displaystyle v_{j}=u_{j}}$ mientras ${displaystyle u_{j}}$ es linealmente independiente de ${displaystyle u_{1},dotscu_{j-1}}$ (y en el caso de que haya tal dependencia entonces uno puede continuar la secuencia escogiendo como ${displaystyle v_{j}}$ un vector arbitrario linealmente independiente ${displaystyle u_{1},dotscu_{j-1}}$ ). Una base que contiene ${displaystyle u_{j}}$ vectores, sin embargo, es probable que estén numéricamente mal condicionados, ya que esta secuencia de vectores es por diseño destinado a converger a un eigenvector de ${displaystyle A}$ . Para evitarlo, se puede combinar la iteración de energía con un proceso Gram-Schmidt, para producir en cambio una base ortonormal de estos subespacios de Krylov.

Elija un vector aleatorio ${displaystyle u_{1}}$ de la norma euclidiana ${displaystyle 1}$ . Vamos. ${displaystyle v_{1}=u_{1}}$ .
Para j=1,... ... ,m− − 1{displaystyle j=1,dotscm-1} ${displaystyle j=1,dotscm-1}$ hacer:
1. Vamos. ${displaystyle u_{j+1}'=Au_{j}}$ .
2. Para todos ${displaystyle k=1,dotscj}$ Deja ${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 de base ${displaystyle v_{1},dotscv_{j}}$ )
3. Vamos. ${displaystyle w_{j+1}=u_{j+1}'-sum _{k=1}^{j}g_{k,j}v_{k}}$ . (Cancelar el componente de ${displaystyle u_{j+1}'}$ que está dentro ${displaystyle operatorname {span} (v_{1},dotscv_{j})}$ )
4. Si ${displaystyle w_{j+1}neq 0}$ Entonces déjalo ${displaystyle u_{j+1}=u_{j+1}'/|u_{j+1}'|}$ y ${displaystyle v_{j+1}=w_{j+1}/|w_{j+1}|}$ ,
  de otro modo ${displaystyle u_{j+1}=v_{j+1}}$ un vector arbitrario de la norma euroclidiana ${displaystyle 1}$ que es ortogonal a todos ${displaystyle v_{1},dotscv_{j}}$ .

La relación entre los vectores de la iteración de poder ${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 no necesitamos realmente ${displaystyle u_{j}}$ vectores para computar estos ${displaystyle v_{j}}$ , porque ${displaystyle u_{j}-v_{j}in operatorname {span} (v_{1},dotscv_{j-1})}$ y por lo tanto la diferencia entre ${displaystyle u_{j+1}'=Au_{j}}$ y ${displaystyle w_{j+1}'=Av_{j}}$ está dentro. ${displaystyle operatorname {span} (v_{1},dotscv_{j})}$ , que es cancelado por el proceso de ortogonalización. Así, la misma base para la cadena de subespacios Krylov es calculada por

Elija un vector aleatorio ${displaystyle v_{1}}$ de la norma euclidiana ${displaystyle 1}$ .
Para j=1,... ... ,m− − 1{displaystyle j=1,dotscm-1} ${displaystyle j=1,dotscm-1}$ hacer:
1. Vamos. ${displaystyle w_{j+1}'=Av_{j}}$ .
2. Para todos ${displaystyle k=1,dotscj}$ Deja ${displaystyle h_{k,j}=v_{k}^{*}w_{j+1}'}$ .
3. Vamos. ${displaystyle w_{j+1}=w_{j+1}'-sum _{k=1}^{j}h_{k,j}v_{k}}$ .
4. Vamos. ${displaystyle h_{j+1,j}=|w_{j+1}|}$ .
5. Si ${displaystyle h_{j+1,j}neq 0}$ Entonces déjalo ${displaystyle v_{j+1}=w_{j+1}/h_{j+1,j}}$ ,
  de otro modo ${displaystyle v_{j+1}}$ un vector arbitrario de la norma euroclidiana ${displaystyle 1}$ que es ortogonal a todos ${displaystyle v_{1},dotscv_{j}}$ .

A priori los coeficientes ${displaystyle h_{k,j}}$ satisfacer satisfacción

{displaystyle Av_{j}=sum _{k=1}^{j+1}h_{k,j}v_{k}}

para todos

<math alttext="{displaystyle jjc)m{displaystyle j) <img alt="{displaystyle j

;

la definición ${displaystyle h_{j+1,j}=|w_{j+1}|}$ puede parecer un poco extraño, pero se ajusta al patrón general ${displaystyle h_{k,j}=v_{k}^{*}w_{j+1}'}$ desde entonces

{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}|.}

Porque el poder iteración vectores ${displaystyle u_{j}}$ que fueron eliminados de esta satisfacción de recursión ${displaystyle u_{j}in operatorname {span} (v_{1},ldotsv_{j}),}$ los vectores ${displaystyle {v_{j}}_{j=1}^{m}}$ y coeficientes ${displaystyle h_{k,j}}$ contener suficiente información ${displaystyle A}$ que todo ${displaystyle u_{1},ldotsu_{m}}$ puede ser calculado, así que nada se perdió cambiando vectores. (De hecho, resulta que los datos reunidos aquí dan unas aproximaciones significativamente mejores del valor más grande que uno obtiene de un número igual de iteraciones en el método de potencia, aunque eso no es necesariamente obvio en este punto.)

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

Elementalmente, ${displaystyle A}$ Hermitian entonces

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

Para $<math alttext="{displaystyle kkc)j- - 1{displaystyle k madej-1} <img alt="{displaystyle k$ sabemos que ${displaystyle Av_{k}in operatorname {span} (v_{1},ldotsv_{j-1})}$ y desde ${displaystyle v_{j}}$ por la construcción es ortogonal a este subespacio, este producto interno debe ser cero. (Esto es esencialmente también la razón por la cual las secuencias de polinomios ortogonales siempre pueden ser dadas una relación de recurrencia de tres plazos.) Para ${displaystyle k=j-1}$ uno se pone

{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 se pone

{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.

Más abstracto, si ${displaystyle V}$ es la matriz con columnas ${displaystyle v_{1},ldotsv_{m}}$ entonces los números ${displaystyle h_{k,j}}$ se puede identificar como elementos de la matriz ${displaystyle H=V^{*}AV}$ , y ${displaystyle h_{k,j}=0}$ para $j+1;}" xmlns="http://www.w3.org/1998/Math/MathML">k■j+1;{displaystyle k títuloj+1;}j+1;}" aria-hidden="true" class="mwe-math-fallback-image-inline mw-invert" src="https://wikimedia.org/api/rest_v1/media/math/render/svg/ac2349a70f6e89d239e891a0836afc1d2c449897" style="vertical-align: -0.671ex; width:9.917ex; height:2.509ex;"/>$ la matriz ${displaystyle H}$ Está arriba de Hessenberg. Desde

{displaystyle H^{*}=left(V^{*}AVright)^{*}=V^{*}A^{*}V=V^{*}AV=H}

la matriz ${displaystyle H}$ es Hermitian. Esto implica que ${displaystyle H}$ es también menor Hessenberg, por lo que debe ser tridiagional. Siendo Hermitian, su diagonal principal es real, y como su primera subdiagonal es real por la construcción, lo mismo es cierto para su primera superdiagonal. Por lo tanto, ${displaystyle H}$ es una matriz real, simétrica— la matriz ${displaystyle T}$ de la especificación del algoritmo Lanczos.

Aproximación simultánea de valores propios extremos

Una manera de caracterizar a los eigenvectores de una matriz hermitiana ${displaystyle A}$ es como puntos estacionarios del cociente Rayleigh

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

En particular, el mayor eigenvalue ${displaystyle lambda _{max }}$ es el máximo global ${displaystyle r}$ y el valor más pequeño ${displaystyle lambda _{min }}$ es el mínimo mundial ${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}$ mínimo ${displaystyle y}$ de ${displaystyle r}$ . Repita eso 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 }$ tales 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}}}

Entonces surge la pregunta de cómo elegir los subespacios para que estas secuencias converjan a un ritmo óptimo.

Desde ${displaystyle x_{j}}$ , la dirección óptima en la que buscar mayores valores ${displaystyle r}$ es el del gradiente ${displaystyle nabla r(x_{j})}$ , y también de ${displaystyle y_{j}}$ la dirección óptima en la que buscar valores más pequeños ${displaystyle r}$ es el 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 lo suficientemente fáciles para calcular en matriz aritmética, pero si uno desea mejorar en ambos ${displaystyle x_{j}}$ y ${displaystyle y_{j}}$ entonces hay dos nuevas direcciones para tener en cuenta: ${displaystyle Ax_{j}}$ y ${displaystyle Ay_{j};}$ desde entonces ${displaystyle x_{j}}$ y ${displaystyle y_{j}}$ pueden ser vectores linealmente independientes (de hecho, están cerca de ortogonal), uno no puede esperar en general ${displaystyle Ax_{j}}$ y ${displaystyle Ay_{j}}$ para 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}}$ son tomados para ser subespacios Krylov, porque entonces ${displaystyle Azin {mathcal {L}}_{j+1}}$ para todos ${displaystyle zin {mathcal {L}}_{j},}$ en particular para ambos ${displaystyle z=x_{j}}$ y ${displaystyle z=y_{j}}$ .

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

{displaystyle {mathcal {L}}_{j}=operatorname {span} (x_{1},Ax_{1},ldotsA^{j-1}x_{1})}

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

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

Desde ${displaystyle j}$ el método de energía iterate ${displaystyle u_{j}}$ pertenece a ${displaystyle {mathcal {L}}_{j},}$ sigue que una iteración para producir ${displaystyle x_{j}}$ y ${displaystyle y_{j}}$ no puede converger más lento que el del método de energía, y logrará más aproximando ambos extremos eigenvalue. Para el subproblema de optimizar ${displaystyle r}$ algunos ${displaystyle {mathcal {L}}_{j}}$ , es conveniente tener una base ortonormal ${displaystyle {v_{1},ldotsv_{j}}}$ para este espacio vectorial. Por lo tanto, estamos nuevamente llevados al problema de calcular iterativamente tal base para la secuencia de subespacios de Krylov.

Convergencia y otras dinámicas

Al analizar la dinámica del algoritmo, es conveniente tomar los eigenvalues y eigenvectores de ${displaystyle A}$ como se da, aunque no sean explícitamente conocidos por el usuario. Para fijar la notación, ${displaystyle lambda _{1}geqslant lambda _{2}geqslant dotsb geqslant lambda _{n}}$ ser los eigenvalues (estos son conocidos por todos ser reales, y por lo tanto posible para ordenar) y dejar ${displaystyle z_{1},dotscz_{n}}$ ser un conjunto ortonormal de eigenvectores tal que ${displaystyle Az_{k}=lambda _{k}z_{k}}$ para todos ${displaystyle k=1,dotscn}$ .

También es conveniente fijar una notación para los coeficientes del vector Lanczos inicial ${displaystyle v_{1}}$ con respecto a esta eigenbasis; ${displaystyle d_{k}=z_{k}^{*}v_{1}}$ para todos ${displaystyle k=1,dotscn}$ Así que ${displaystyle textstyle v_{1}=sum _{k=1}^{n}d_{k}z_{k}}$ . Un vector inicial ${displaystyle v_{1}}$ El agotamiento de algunos eigencomponentes retrasará la convergencia al eigenvalue correspondiente, e incluso aunque esto simplemente salga como un factor constante en los límites del error, el agotamiento sigue siendo indeseable. Una técnica común para evitar ser golpeada constantemente por él es elegir ${displaystyle v_{1}}$ por primera vez dibujar los elementos aleatoriamente según la misma distribución normal con media ${displaystyle 0}$ y luego reescala el vector a la norma ${displaystyle 1}$ . Antes de la escalada, esto causa los coeficientes ${displaystyle d_{k}}$ ser independiente normalmente distribuida variables estocásticas de la misma distribución normal (ya que el cambio de coordenadas es unitario), y después de escalar el vector ${displaystyle (d_{1},dotscd_{n})}$ tendrá una distribución uniforme en el ámbito de la unidad en ${displaystyle mathbb {C} ^{n}}$ . Esto hace posible ligar la probabilidad de que por ejemplo $<math alttext="{displaystyle |d_{1}|Silenciod1Silencioc)ε ε {displaystyle Silencio.<img alt="{displaystyle |d_{1}|$ .

El hecho de que el algoritmo de Lanczos es coordinámica – las operaciones sólo miran los productos interiores de vectores, nunca en elementos individuales de vectores – hace fácil construir ejemplos con eigenestructura conocida para ejecutar el algoritmo en: ${displaystyle A}$ una matriz diagonal con los eigenvalues deseados en la diagonal; siempre y cuando el vector inicial ${displaystyle v_{1}}$ tiene suficientes elementos no cero, el algoritmo producirá una matriz simétrica tridiagonal general como ${displaystyle T}$ .

Teoría de la convergencia de Kaniel-Paige

Después ${displaystyle m}$ pasos de la iteración del algoritmo Lanczos, ${displaystyle T}$ es un ${displaystyle mtimes m}$ matriz simétrica real, que similarmente a lo anterior tiene ${displaystyle m}$ eigenvalues ${displaystyle theta _{1}geqslant theta _{2}geqslant dots geqslant theta _{m}.}$ Por convergencia se entiende principalmente la convergencia ${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 en segundo lugar la convergencia de algún rango ${displaystyle theta _{1},ldotstheta _{k}}$ of eigenvalues of ${displaystyle T}$ a sus homólogos ${displaystyle lambda _{1},ldotslambda _{k}}$ de ${displaystyle A}$ . La convergencia para el algoritmo Lanczos es a menudo órdenes de magnitud más rápido que eso para el algoritmo de iteración de energía.

Los límites para ${displaystyle theta _{1}}$ proviene de la interpretación anterior de eigenvalues como valores extremos del cociente Rayleigh ${displaystyle r(x)}$ . Desde ${displaystyle lambda _{1}}$ es a priori el máximo ${displaystyle r}$ en todo el ${displaystyle mathbb {C} ^{n},}$ mientras que ${displaystyle theta _{1}}$ es simplemente el máximo en un ${displaystyle m}$ -dimensional Subespacial de Krylov, trivialmente obtenemos ${displaystyle lambda _{1}geqslant theta _{1}}$ . Por el contrario, cualquier punto ${displaystyle x}$ en que el subespacio Krylov proporciona un límite inferior ${displaystyle r(x)}$ para ${displaystyle theta _{1}}$ , así que si se puede exhibir un punto ${displaystyle lambda _{1}-r(x)}$ es pequeño entonces esto proporciona un límite apretado ${displaystyle theta _{1}}$ .

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

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

para que cualquier elemento de ella pueda expresarse como ${displaystyle p(A)v_{1}}$ para algunos polinomios ${displaystyle p}$ de grado en la mayoría ${displaystyle m-1}$ ; los coeficientes de ese polinomio son simplemente los coeficientes en la combinación lineal de los vectores ${displaystyle v_{1},Av_{1},A^{2}v_{1},ldotsA^{m-1}v_{1}}$ . El polinomio que queremos resultará tener coeficientes reales, pero por el momento debemos permitir también los coeficientes complejos, y escribiremos ${displaystyle p^{*}}$ para el polinomio obtenido por complejo conjugando todos los coeficientes de ${displaystyle p}$ . En esta parametrización del subespacio 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 eigenvectores, obtenemos

{displaystyle Av_{1}=Asum _{k=1}^{n}d_{k}z_{k}=sum _{k=1}^{n}d_{k}lambda _{k}z_{k}}

y más en general

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

para cualquier polinomio ${displaystyle q}$ .

Así

{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 el ${displaystyle k=1}$ el término desaparece en el numerador, pero no en el denominador. Así si uno puede elegir ${displaystyle p}$ para ser grande ${displaystyle lambda _{1}}$ pero pequeño en todos los otros eigenvalues, uno conseguirá un vínculo estrecho en el error ${displaystyle lambda _{1}-theta _{1}}$ .

Desde ${displaystyle A}$ tiene muchos más valores que ${displaystyle p}$ tiene coeficientes, esto puede parecer un orden alto, pero una forma de conocerlo es utilizar polinomios Chebyshev. Escritura ${displaystyle c_{k}}$ para el grado ${displaystyle k}$ Polinomio Chebyshev del primer tipo (lo que satisfice ${displaystyle c_{k}(cos x)=cos(kx)}$ para todos ${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 ella. Con alguna escalada del argumento, podemos tener que mapear todos los eigenvalues excepto ${displaystyle lambda _{1}}$ en ${displaystyle [-1,1]}$ . Vamos.

{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}}$ , utilizar en su lugar el mayor eigenvalue estrictamente menos que ${displaystyle lambda _{1}}$ ), entonces el valor máximo de ${displaystyle |p(lambda _{k})|^{2}}$ para ${displaystyle kgeqslant 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 {2lambda _{1}-lambda _{2}-lambda _{n}}{lambda _{2}-lambda _{n}}}right)=c_{m-1}left(2{frac {lambda _{1}-lambda _{2}}{lambda _{2}-lambda _{n}}}+1right);}

la cantidad

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

(es decir, la relación entre el primer espacio propio y el diámetro del resto del espectro) es, por tanto, de importancia clave para la tasa de convergencia. tambien escribiendo

{displaystyle R=e^{operatorname {arcosh} (1+2rho)}=1+2rho +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}(2rho +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+2rho))}}\[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}}$ para cada iteración extra.

Para la comparación, se puede considerar cómo la tasa de convergencia del método de potencia depende de ${displaystyle rho }$ , pero como el método de poder es sensible principalmente al cociente entre valores absolutos de los eigenvalues, necesitamos ${displaystyle |lambda _{n}|leqslant |lambda _{2}|}$ para el eigengap entre ${displaystyle lambda _{1}}$ y ${displaystyle lambda _{2}}$ ser el dominante. Bajo esa limitación, el caso que más favorece el método de potencia es que ${displaystyle lambda _{n}=-lambda _{2}}$ , así que considera eso. Tarde en el 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},}

donde cada nueva iteración multiplica eficazmente ${displaystyle z_{2}}$ -La libertad ${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+2rho }}.}

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 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+2rho)^{-2}}$ para cada iteración. La diferencia así se reduce a eso entre ${displaystyle 1+2rho }$ y ${displaystyle R=1+2rho +2{sqrt {rho ^{2}+rho }}}$ . En el ${displaystyle rho gg 1}$ región, este último es más como ${displaystyle 1+4rho }$ , y realiza como el método de energía con un eigengap dos veces más grande; una mejora notable. El caso más difícil es sin embargo el de ${displaystyle rho ll 1,}$ en que ${displaystyle Rapprox 1+2{sqrt {rho }}}$ es una mejora aún mayor en el eigengap; el ${displaystyle rho gg 1}$ región es donde el algoritmo de la convergencia de Lanczos hace que más pequeña mejora en el método de potencia.

Estabilidad numérica

La estabilidad significa 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.

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

Los usuarios de este algoritmo deben poder encontrar y eliminar esos mensajes "falsos" valores propios. Las implementaciones prácticas del algoritmo Lanczos van en tres direcciones para combatir este problema de estabilidad:

Prevenir la pérdida de ortogonalidad,
Recuperar la ortogonalidad después de que se genere la base.
Después de que se identifiquen los eigenvalues buenos y "espuriosos", retire los espurios.

Variaciones

Existen variaciones en el algoritmo de Lanczos donde los vectores involucrados son matrices altas y estrechas en lugar de vectores y las constantes de normalización son matrices cuadradas pequeñas. Estos se denominan "bloque" Los algoritmos de Lanczos 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 un cierto número de iteraciones. Una de las variaciones reiniciadas más influyentes es el método Lanczos reiniciado implícitamente, que se implementa en ARPACK. Esto ha llevado a una serie de otras variaciones reiniciadas, como la bidiagonalización reiniciada de Lanczos. Otra variación reiniciada exitosa es el método Lanczos de reinicio grueso, que se implementó en un paquete de software llamado TRLan.

Espacio nulo sobre un campo 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 grandes matrices dispersas sobre campos finitos y el conjunto de personas interesadas en grandes problemas de valores propios apenas se superponen, esto a menudo también se denomina algoritmo de bloque de Lanczos sin causar una confusión irrazonable.

Aplicaciones

Los algoritmos de 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 de largo plazo implementan esta operación, el algoritmo Lanczos se puede aplicar eficientemente a los documentos de texto (ver indexación semántica latente). Eigenvectors también son importantes para métodos de clasificación a gran escala como el algoritmo HITS desarrollado por Jon Kleinberg, o el algoritmo PageRank utilizado por Google.

Los algoritmos de Lanczos también se utilizan en física de la materia condensada como método para resolver hamiltonianos de sistemas de electrones fuertemente correlacionados, así como en códigos de modelos de capas en física nuclear.

Implementaciones

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

MATLAB y GNU Octave vienen con ARPACK integrado. 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.

Hay disponible una implementación en Matlab del algoritmo Lanczos (tenga en cuenta los problemas de precisión) como parte del paquete Matlab de propagación de creencias gaussianas. La biblioteca de filtrado colaborativo GraphLab incorpora una implementación paralela a gran escala del algoritmo Lanczos (en C++) para multinúcleo.

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

Contenido relacionado

Más resultados...