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 método iterativo . Arnoldi encuentra una aproximación a los valores propios y vectores propios de matrices generales (posiblemente no hermitianas ) mediante la construcción de una base ortonormal del subespacio de Krylov , lo que lo hace particularmente útil cuando se trata de 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, a diferencia de los llamados métodos directos que deben completarse para dar resultados útiles (ver, por ejemplo, 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 hermitianas 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 poder.

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, muchos cálculos potencialmente útiles se desperdician al utilizar sólo el resultado final . Esto sugiere que, en cambio, formamos la llamada 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 resultante de vectores 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 se aproxima al 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 , tal que para cada n , los vectores q ₁ ,... , q _n abarca 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 j -loop proyecta el componente de en las direcciones de . Esto asegura la ortogonalidad de todos los vectores generados. ${\displaystyle q_{k}}$ ${\displaystyle q_{1},\dots ,q_{k-1}}$

El algoritmo falla 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 valores propios siguiente 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 ): """Calcular 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:  vector inicial similar a una matriz (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 (n + 1) xn. A sobre la base Q. Es el Alto Hessenberg.  """ 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 . norma ( b , 2 ) # Úselo 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 en 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 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.  devolver  Q ,  h  devolver  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 (de Hessenberg superior ) formada por los números h _{j , k} calculado por el algoritmo:

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

El método de ortogonalización debe elegirse específicamente de manera que los componentes inferiores de Arnoldi/Krylov se eliminen de los vectores superiores de Krylov. 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 se puede caracterizar por la siguiente condición de optimización. El polinomio característico de H _n minimiza || pag ( A ) q ₁ || ₂ entre todos los polinomios mónicos de grado n . Este problema de optimización tiene una solución única si y sólo si la iteración de Arnoldi no se descompone. ${\displaystyle {\mathcal {K}}_{n}}$ ${\displaystyle {\mathcal {K}}_{n}}$

La relación entre las matrices Q en iteraciones posteriores está 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 agregando una fila adicional a H _n .

Encontrar valores propios con la iteración de Arnoldi

La idea de la iteración de Arnoldi como 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 algo relacionado, el algoritmo de Francis . Además, se puede considerar que el algoritmo de Francis está relacionado con iteraciones de potencia, que operan en el subespacio anidado de Krylov. De hecho, la forma más básica del algoritmo de Francis parece ser elegir b para que sea igual a Ae ₁ y extender n a la dimensión completa de A. Las versiones mejoradas incluyen uno o más turnos y se pueden aplicar potencias más altas de A en un solo paso. ^[2]

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

A menudo se observa en la práctica que algunos de los valores propios de Ritz convergen a valores propios de A. Dado que H _n es n -by- n , tiene como máximo n valores propios y no todos los valores propios de A pueden aproximarse. Normalmente, los valores propios de Ritz convergen a los valores propios más grandes de A. Para obtener los valores propios más pequeños de A , se debe utilizar la (operación) inversa de A. Esto puede estar relacionado con la caracterización de H _n como la matriz cuyo polinomio característico minimiza || pag ( A ) q ₁ || de la siguiente manera. Una buena manera 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 valor propio de A. Por lo tanto, los ceros de p (y por tanto los valores propios de Ritz) estarán cerca de los valores propios de A.

Sin embargo, los detalles aún no se comprenden completamente. Esto contrasta con el caso en el que A es hermitiano . 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) con los valores propios (negro) de una matriz de 400x400, compuesta de valores aleatorios uniformes en el dominio [-0,5 +0,5]

Iteración Arnoldi reiniciada

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

Ver 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, NOSOTROS (1951). "El principio de iteraciones minimizadas en la solución del problema de valores propios de matrices". Trimestral de Matemática Aplicada . 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 y Aplicaciones de Matrices . 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 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 y Aplicaciones de Matrices . 23 (3): 601–614. doi :10.1137/S0895479800371529. ISSN  0895-4798.

• WE Arnoldi, "El principio de iteraciones minimizadas en la solución del problema de valores propios de la matriz", Quarterly of Applied Mathematics , volumen 9, páginas 17-29, 1951.
• Yousef Saad , Métodos numéricos para grandes problemas de 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 precondicionados 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().