Resolver un sistema con una actualización diagonal de rango pequeño

Supongamos que tengo el sistema lineal original grande y disperso: . Ahora, no tengo A - 1 ya que A es demasiado grande para factorizar o cualquier tipo de descomposición de A , pero supongo que tengo la solución x 0 encontrada con una solución iterativa.$A\textbf{x}_0=\textbf{b}_0$$A^{-1}$$A$$\textbf{x}_0$

Ahora, deseo aplicar una actualización de rango pequeño a la diagonal de A (cambiar algunas de las entradas diagonales): donde D es una matriz diagonal con mayormente 0 en su diagonal y algunas valores distintos de cero. Si tuviera A - 1 podría aprovechar la fórmula de Woodbury para aplicar una actualización a la inversa. Sin embargo, no tengo esto disponible. ¿Hay algo que pueda hacer aparte de resolver todo el sistema nuevamente? ¿Hay alguna forma de que pueda idear un preacondicionador M que sea fácil \ más fácil de invertir, de modo que M A 1$(A+D)\textbf{x}_1=\textbf{b}_0$$D$$A^{-1}$$M$ , de modo que todo lo que tendría que hacer si tengo x 0 es aplicar M - 1 y un método iterativo convergería en un par / pocas iteraciones.$MA_1 \approx A_0$$\textbf{x}_0$$M^{-1}$

1. Guardar en las columnas de dos matrices y C todos los vectores$B$$C$ a los que aplicó la matriz en las iteraciones anteriores y los resultados c j = A b j .$b_j$$c_j=Ab_j$

2. Para cada nuevo sistema (o A x = b ′ , que es el caso especial D = 0 ), resuelva aproximadamente el sistema lineal sobredeterminado ( C + D B ) y ≈ b ′ , p. Ej. , seleccionando un subconjunto de las filas (posiblemente todas) y utilizando un método denso mínimo cuadrado. Tenga en cuenta que solo la parte seleccionada de C + D B necesita ser ensamblada; ¡Esta es una operación rápida!$(A+D)x'=b'$$Ax=b'$$D=0$$(C+DB)y\approx b'$$C+DB$

3. Pon . Esta es una buena aproximación inicial con la cual comenzar la iteración para resolver ( A + D ) $x_0=By$ . En caso de que se deban procesar más sistemas, use los productos de vector de matriz en esta nueva iteración para extender las matrices B y C en el subsistema resultante.$(A+D)x'=b'$$B$$C$

Si las matrices y C no caben en la memoria principal, almacene B en el disco y seleccione el subconjunto de filas por adelantado. Esto le permite mantener en el núcleo la parte relevante de B y C necesaria para formar el sistema de mínimos cuadrados, y el siguiente x 0 puede calcularse mediante una pasada a través de B con poco uso de la memoria central.$B$$C$$B$$B$$C$$x_0$$B$

Las filas deben seleccionarse de tal manera que correspondan aproximadamente a una discreta discreción del problema completo. Tomar cinco veces más filas que el número total de multiplicaciones esperadas de vectores de matriz debería ser suficiente.

Editar: ¿Por qué funciona esto? Por construcción, las matrices y C están relacionados por C = A B . Si el subespacio atravesado por las columnas de B contiene el vector de solución exacto x ' (una situación rara pero simple), entonces x ' tiene la forma x ' = B y para algunos y . Sustituyendo esto en la ecuación que define x ′ se obtiene la ecuación ( C + D B ) y = b $B$$C$$C=AB$$B$$x'$$x'$$x'=By$$y$$x'$$(C+DB)y= b'$. Por lo tanto, en este caso, el proceso anterior da como punto de partida , que es la solución exacta.$x_0=By=x'$

En general, no se puede esperar encuentre en el espacio de la columna de B , pero el punto de partida generado será el punto en este espacio nublado más cercano a x ' , en una métrica determinada por las filas seleccionadas. Por lo tanto, es probable que sea una aproximación sensata. A medida que se procesan más sistemas, el espacio de la columna crece y la aproximación probablemente mejorará mucho, por lo que uno puede esperar converger en cada vez menos iteraciones.$x'$$B$$x'$

Edit2: Acerca del subespacio generado: si uno resuelve cada sistema con un método de Krylov, los vectores utilizados para obtener el punto de partida para el segundo sistema abarcan el subespacio de Krylov del primer lado derecho. Por lo tanto, se obtiene una buena aproximación cada vez que este subespacio de Krylov contiene un vector cercano a la solución de su segundo sistema. En general, los vectores utilizados para obtener el punto de partida para el st abarcan un espacio que contiene el subespacio de Krylov de los primeros k lados derechos.$(k+1)$$k$