Iteración de Arnoldi

Método iterativo para aproximar vectores propios

En álgebra lineal numérica , la iteración de Arnoldi es un algoritmo de valores propios y un ejemplo importante de un método iterativo . Arnoldi encuentra una aproximación a los valores propios y vectores propios de matrices generales (posiblemente no hermíticas ) mediante la construcción de una base ortonormal del subespacio de Krylov , lo que lo hace particularmente útil cuando se trabaja con matrices dispersas grandes .

El método de Arnoldi pertenece a una clase de algoritmos de álgebra lineal que dan un resultado parcial después de un pequeño número de iteraciones, en contraste con los llamados métodos directos que deben completarse para dar resultados útiles (véase por ejemplo la transformación de Householder ). El resultado parcial en este caso son los primeros vectores de la base que está construyendo el algoritmo.

Cuando se aplica a matrices hermíticas, se reduce al algoritmo de Lanczos . La iteración de Arnoldi fue inventada por WE Arnoldi en 1951. ^[1]

Subespacios de Krylov y la iteración de potencia

Un método intuitivo para encontrar el valor propio más grande (en valor absoluto) de una matriz m × m dada es la iteración de potencia : comenzando con un vector inicial arbitrario b , calcule Ab , A ² b , A ³ b , ... normalizando el resultado después de cada aplicación de la matriz A . ${\displaystyle A}$

Esta secuencia converge al vector propio correspondiente al valor propio con el valor absoluto más grande, . Sin embargo, se desperdician muchos cálculos potencialmente útiles si se utiliza únicamente el resultado final, . Esto sugiere que, en su lugar, formamos la denominada matriz de Krylov : ${\displaystyle \lambda _{1}}$ ${\displaystyle A^{n-1}b}$

${\displaystyle K_{n}={\begin{bmatrix}b&Ab&A^{2}b&\cdots &A^{n-1}b\end{bmatrix}}.}$

Las columnas de esta matriz no son en general ortogonales , pero podemos extraer una base ortogonal , mediante un método como la ortogonalización de Gram-Schmidt . El conjunto de vectores resultante es, por tanto, una base ortogonal del subespacio de Krylov , . Podemos esperar que los vectores de esta base abarquen buenas aproximaciones de los vectores propios correspondientes a los valores propios más grandes, por la misma razón que aproxima el vector propio dominante. ${\displaystyle {\mathcal {K}}_{n}}$ ${\displaystyle n}$ ${\displaystyle A^{n-1}b}$

La iteración de Arnoldi

La iteración de Arnoldi utiliza el proceso de Gram-Schmidt modificado para producir una secuencia de vectores ortonormales, q ₁ , q ₂ , q ₃ , ..., llamados vectores de Arnoldi , de modo que para cada n , los vectores q ₁ , ..., q _n abarcan el subespacio de Krylov . Explícitamente, el algoritmo es el siguiente: ${\displaystyle {\mathcal {K}}_{n}}$

Comience con un vector arbitrario q 1 con norma 1. Repita para  k = 2, 3, ... q k  := A  q k −1  para  j  de 1 a  k − 1 h j , k −1  := q j *  q k  q k  := q k − h j , k −1  q j  h k , k −1  := ‖ q k ‖ q k  := q k / h k , k −1

El bucle j proyecta el componente de en las direcciones de . Esto garantiza la ortogonalidad de todos los vectores generados. ${\displaystyle q_{k}}$ ${\displaystyle q_{1},\dots ,q_{k-1}}$

El algoritmo deja de funcionar cuando q _k es el vector cero. Esto sucede cuando el polinomio mínimo de A es de grado k . En la mayoría de las aplicaciones de la iteración de Arnoldi, incluido el algoritmo de valor propio que se muestra a continuación y GMRES , el algoritmo ha convergido en este punto.

Cada paso del bucle k requiere un producto matriz-vector y aproximadamente 4 mk operaciones de punto flotante.

En el lenguaje de programación Python con soporte de la biblioteca NumPy :

importar  numpy  como  npdef  arnoldi_iteration ( A ,  b ,  n :  int ): """Calcula una base del subespacio (n + 1)-Krylov de la matriz A.  Este es el espacio abarcado por los vectores {b, Ab, ..., A^nb}. Parámetros  ----------  A : array_like  Una matriz m × m.  b : array_like  Vector inicial (longitud m).  n : int  Uno menos que la dimensión del subespacio de Krylov, o equivalentemente el *grado* del espacio de Krylov. Debe ser >= 1.  Devuelve  -------  Q : numpy.array  Una matriz mx (n + 1), donde las columnas son una base ortonormal del subespacio de Krylov.  h : numpy.array  Una matriz xn (n + 1). Una sobre base Q. Es Hessenberg superior.  """ eps = 1e-12 h = np . zeros (( n + 1 , n )) Q = np . zeros (( A . shape [ 0 ], n + 1 )) # Normalizar el vector de entrada Q [:, 0 ] = b / np . linalg . norm ( b , 2 ) # Usarlo como el primer vector de Krylov para k en el rango ( 1 , n + 1 ): v = np . dot ( A , Q [:, k - 1 ]) # Generar un nuevo vector candidato para j en el rango ( k ): # Restar las proyecciones sobre los vectores anteriores h [ j , k - 1 ] = np . dot ( Q [:, j ] . conj (), v ) v = v - h [ j , k - 1 ] * Q [:, j ] h [ k , k - 1 ] = np . linalg . norma ( v , 2 ) si h [ k , k - 1                                                                            ]  >  eps :  # Agrega el vector producido a la lista, a menos que  Q [:,  k ]  =  v  /  h [ k ,  k  -  1 ]  else :  # Si eso sucede, deja de iterar.  return  Q ,  h  return  Q ,  h

Propiedades de la iteración de Arnoldi

Sea Q _n la matriz m por n formada por los primeros n vectores de Arnoldi q ₁ , q ₂ , ..., q _n , y sea H _n la matriz ( Hessenberg superior ) formada por los números h _{j , k} calculados por el algoritmo:

${\displaystyle H_{n}=Q_{n}^{*}AQ_{n}.}$

El método de ortogonalización debe elegirse específicamente de modo que los componentes de Arnoldi/Krylov inferiores se eliminen de los vectores de Krylov superiores. Como se puede expresar en términos de q ₁ , ..., q _{i +1} por construcción, son ortogonales a q _{i +2} , ..., q _n , ${\displaystyle Aq_{i}}$

Entonces tenemos

${\displaystyle H_{n}={\begin{bmatrix}h_{1,1}&h_{1,2}&h_{1,3}&\cdots &h_{1,n}\\h_{2,1}&h_{2,2}&h_{2,3}&\cdots &h_{2,n}\\0&h_{3,2}&h_{3,3}&\cdots &h_{3,n}\\\vdots &\ddots &\ddots &\ddots &\vdots \\0&\cdots &0&h_{n,n-1}&h_{n,n}\end{bmatrix}}.}$

La matriz H _n puede verse como A en el subespacio con los vectores de Arnoldi como base ortogonal; A se proyecta ortogonalmente sobre . La matriz H _n puede caracterizarse por la siguiente condición de optimalidad. El polinomio característico de H _n minimiza || p ( A ) q ₁ || ₂ entre todos los polinomios mónicos de grado n . Este problema de optimalidad tiene una solución única si y solo si la iteración de Arnoldi no se rompe. ${\displaystyle {\mathcal {K}}_{n}}$ ${\displaystyle {\mathcal {K}}_{n}}$

La relación entre las matrices Q en iteraciones posteriores viene dada por

${\displaystyle AQ_{n}=Q_{n+1}{\tilde {H}}_{n}}$

dónde

${\displaystyle {\tilde {H}}_{n}={\begin{bmatrix}h_{1,1}&h_{1,2}&h_{1,3}&\cdots &h_{1,n}\\h_{2,1}&h_{2,2}&h_{2,3}&\cdots &h_{2,n}\\0&h_{3,2}&h_{3,3}&\cdots &h_{3,n}\\\vdots &\ddots &\ddots &\ddots &\vdots \\\vdots &&0&h_{n,n-1}&h_{n,n}\\0&\cdots &\cdots &0&h_{n+1,n}\end{bmatrix}}}$

es una matriz ( n +1) por n formada añadiendo una fila extra a H _n .

Encontrar valores propios con la iteración de Arnoldi

La idea de la iteración de Arnoldi como un algoritmo de valores propios es calcular los valores propios en el subespacio de Krylov. Los valores propios de H _n se denominan valores propios de Ritz . Dado que H _n es una matriz de Hessenberg de tamaño modesto, sus valores propios se pueden calcular de manera eficiente, por ejemplo, con el algoritmo QR o, de alguna manera relacionado, el algoritmo de Francis . También se puede considerar que el algoritmo de Francis en sí está relacionado con las iteraciones de potencia, que operan en el subespacio de Krylov anidado. De hecho, la forma más básica del algoritmo de Francis parece ser elegir b como igual a Ae ₁ y extender n a la dimensión completa de A. Las versiones mejoradas incluyen uno o más desplazamientos, y se pueden aplicar potencias mayores de A en un solo paso. ^[2]

Este es un ejemplo del método Rayleigh-Ritz .

En la práctica, se observa a menudo que algunos de los autovalores de Ritz convergen a los autovalores de A. Dado que H _n es n -por- n , tiene como máximo n autovalores, y no todos los autovalores de A pueden aproximarse. Normalmente, los autovalores de Ritz convergen a los mayores autovalores de A. Para obtener los menores autovalores de A , se debe utilizar en su lugar la operación inversa de A. Esto se puede relacionar con la caracterización de H _n como la matriz cuyo polinomio característico minimiza || p ( A ) q ₁ || de la siguiente manera. Una buena forma de hacer que p ( A ) sea pequeño es elegir el polinomio p tal que p ( x ) sea pequeño siempre que x sea un autovalor de A. Por lo tanto, los ceros de p (y, por tanto , los autovalores de Ritz) estarán cerca de los autovalores de A.

Sin embargo, los detalles aún no se comprenden por completo. Esto contrasta con el caso en el que A es hermítico . En esa situación, la iteración de Arnoldi se convierte en la iteración de Lanczos , para la cual la teoría es más completa.

Iteración de Arnoldi que demuestra la convergencia de los valores de Ritz (rojo) a los valores propios (negro) de una matriz de 400x400, compuesta de valores aleatorios uniformes en el dominio [-0,5 +0,5]

Se reinició la iteración de Arnoldi

Debido a consideraciones prácticas de almacenamiento, las implementaciones comunes de los métodos de Arnoldi suelen reiniciarse después de un número fijo de iteraciones. Un enfoque es el método de Arnoldi reiniciado implícitamente (IRAM) ^[3] de Lehoucq y Sorensen, que se popularizó en el paquete de software libre y de código abierto ARPACK . ^[4] Otro enfoque es el algoritmo Krylov-Schur de GW Stewart, que es más estable y más simple de implementar que el IRAM. ^[5]

Véase también

El método residual mínimo generalizado (GMRES) es un método para resolver Ax = b basado en la iteración de Arnoldi.

Referencias

1. Arnoldi, WE (1951). "El principio de iteraciones minimizadas en la solución del problema de valores propios de matrices". Quarterly of Applied Mathematics . 9 (1): 17–29. doi : 10.1090/qam/42792 . ISSN  0033-569X.
2. David S. Watkins. Algoritmo de Francis, Universidad Estatal de Washington. Consultado el 14 de diciembre de 2022
3. RB Lehoucq y DC Sorensen (1996). "Técnicas de deflación para una iteración de Arnoldi reiniciada implícitamente". Revista SIAM sobre análisis de matrices y aplicaciones . 17 (4): 789–821. doi :10.1137/S0895479895281484. hdl : 1911/101832 .
4. RB Lehoucq; DC Sorensen y 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. Archivado desde el original el 26 de junio de 2007. Consultado el 30 de junio de 2007 .
5. Stewart, GW (2002). "Un algoritmo de Krylov--Schur para grandes problemas propios". Revista SIAM sobre análisis de matrices y aplicaciones . 23 (3): 601–614. doi :10.1137/S0895479800371529. ISSN  0895-4798.

• WE Arnoldi, "El principio de iteraciones minimizadas en la solución del problema del valor propio de la matriz", Quarterly of Applied Mathematics , volumen 9, páginas 17–29, 1951.
• Yousef Saad , Métodos numéricos para problemas de grandes valores propios , Manchester University Press, 1992. ISBN 0-7190-3386-1 . 
• Lloyd N. Trefethen y David Bau, III, Álgebra lineal numérica , Sociedad de Matemáticas Industriales y Aplicadas, 1997. ISBN 0-89871-361-7 . 
• Jaschke, Leonhard: Métodos de Arnoldi preacondicionados para sistemas de ecuaciones no lineales . (2004). ISBN 2-84976-001-3 
• Implementación: Matlab viene con ARPACK incorporado. Tanto las matrices almacenadas como las implícitas se pueden analizar mediante la función eigs().