Iteración de Arnoldi

En álgebra lineal numérica , la iteración de Arnoldi es un algoritmo para calcular valores propios . Arnoldi encuentra una aproximación de los valores propios y vectores propios de matrices generales (posiblemente no hermitianas ) mediante la construcción de una base ortonormal del subespacio de Krylov .

El método de Arnoldi se refiere a algoritmos de álgebra lineal que permiten obtener una solución parcial después de un pequeño número de iteraciones, a diferencia de los llamados métodos directos, que deben completarse por completo para producir resultados satisfactorios (por ejemplo , la reflexión del jefe de hogar ).

Si el algoritmo se aplica sobre matrices hermitianas, entonces se reduce al algoritmo de Lanczos . La iteración de Arnoldi fue inventada por Walter Edwin Arnoldi en 1951.

El subespacio de Krylov y el método del poder

Un método intuitivo para encontrar el valor propio más grande (módulo) de una matriz dada de dimensiones es el método de la potencia : comience con un vector inicial arbitrario , calcule y normalice el resultado después de cada cálculo. $A$ $m \veces m$ $b$ $Ab,A^{2}b,A^{3}b,\puntos$

Esta secuencia converge al vector propio del valor propio correspondiente con módulo máximo. Esto sugiere que se desperdicia una gran cantidad de cálculo, ya que al final, solo se usa el resultado final . Entonces, en cambio, compongamos la llamada matriz de Krylov : $\lambda _{1}$ ${\ estilo de visualización A ^ {n-1} b}$

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

Las columnas de esta matriz generalmente no son ortogonales, pero podemos obtener una base ortogonal de ellas utilizando la ortogonalización de Gram-Schmidt . El conjunto de vectores resultante será la base ortogonal del subespacio de Krylov . Se puede esperar que los vectores de esta base sean una buena aproximación a los vectores correspondientes a los valores propios más grandes en valor absoluto. ${\mathcal {K}}_{n}$ $norte$

Iteración de Arnoldi

La iteración de Arnoldi usa un proceso de Gram-Schmidt estabilizado para obtener una secuencia de vectores ortonormales , llamados vectores de Arnoldi , tales que, para cada uno, los vectores son la base de un subespacio de Krylov . El algoritmo se ve así: ${\ estilo de visualización q_ {1}, q_ {2}, q_ {3}, \ puntos}$ $norte$ ${\displaystyle q_{1},\puntos,q_{n))$ ${\mathcal {K}}_{n}$

Comience con un vector arbitrario q 1 con norma 1.
Repita para k = 2, 3, …
- $q_{k}\leftarrow Aq_{k-1}\,$
- para j = 1 ... k − 1
  - $h_{j,k-1}\leftarrow q_{j}^{*}q_{k}\,$
  - $q_{k}\leftarrow q_{k}-h_{j,k-1}q_{j}\,$
- $h_{k,k-1}\flecha izquierda \|q_{k}\|\,$
- $q_{k}\leftarrow {\frac {q_{k}}{h_{k,k-1}}}\,$

El loop over proyecta el componente sobre Esto asegura la ortogonalidad de todos los vectores construidos. $j$ $q_{k}$ $q_{1},\puntos,q_{k-1}.$

El algoritmo se detiene cuando es un vector cero. Esto sucede cuando el polinomio mínimo de la matriz será de grado $q_{k}$ $A$ $k.$

Cada paso del ciclo for realiza una multiplicación matriz-vector y operaciones fraccionarias. $k$ ${\ estilo de visualización 4mk}$

En lenguaje de programación python:

importar numpy como np def arnoldi_iteration ( A , b , n : int ): """Calcula una base del (n + 1)-subespacio de Krylov de A: el espacio abarcado por {b, Ab, ..., A^nb}. Argumentos A: matriz m × m b: vector inicial (longitud m) n: dimensión del subespacio de Krylov, debe ser >= 1 Devuelve la matriz Q: mx (n + 1), las columnas son una base ortonormal del subespacio de Krylov. h: (n + 1) matriz xn, A sobre base Q. Es Hessenberg superior. """ m = A . forma [ 0 ] h = np . ceros (( n + 1 , n )) Q = np . ceros (( m , n + 1 )) q = b / np . linalg . norma ( b ) # Normalizar el vector de entrada Q [:, 0 ] = q # Usarlo como el primer vector de Krylov para k en el rango ( n ): v = A . punto ( q ) # Generar un nuevo vector candidato para j en el rango ( k + 1 ): # Restar las proyecciones de los vectores anteriores h [ j , k ] = np . punto ( Q [:, j ] . conj (), v ) v = v - h [ j , k ] * Q [:, j ] h [ k + 1 , k ] = np . Linalg . norma ( v ) eps = 1e-12 # Si v es más corto que este umbral, es el vector cero si h [ k + 1 , k ] > eps : # Agregue el vector producido a la lista, a menos que q = v / h [ k + 1 , k ] # se produce el vector cero. Q [:, k + 1 ] = q else : # Si eso sucede, deja de iterar. devolver Q , h devolver Q , h