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

9

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.UNAX0 0=si0 0UNA-1UNAX0 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(UNA+re)X1=si0 0reUNA-1METRO , 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.METROUNA1UNA0 0X0 0METRO-1

Costis
fuente
¿Estás comenzando con un buen preacondicionador para y quieres saber cómo actualizarlo? ¿Qué rango tiene la actualización? (Una actualización de rango 1000 es "pequeña" en comparación con una matriz de tamaño 10 9 pero no pequeña en términos de recuento de iteraciones.)UNA1000109 9
Jed Brown
tiene un tamaño aproximado de 10 6 a 10 7 , y la actualización es <1000 (probablemente <100) elementos. Estoy usando un tipo de preacondicionador diagonal para A que funciona muy bien, por lo que actualizarlo sería trivial, pero me preguntaba si hay algo mejor que pueda hacer en lugar de resolver el nuevo sistema desde cero. UNA106107 7
Costis
2
La solución de un sistema no le dice mucho al respecto. Si resuelve el mismo sistema varias veces, el mapa inverso en esos vectores (y / o espacios de Krylov asociados) le brinda información que puede usarse para acelerar la convergencia. ¿Cuántos sistemas estás resolviendo en cada caso?
Jed Brown
Actualmente sólo soy la solución para una RHS ( vectorial) con cada una matriz antes de modificar una . siUNAUNA
Costis

Respuestas:

4
  1. Guardar en las columnas de dos matrices y C todos los vectoressiC a los que aplicó la matriz en las iteraciones anteriores y los resultados c j = A b j .sijCj=UNAsij

  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!(UNA+re)X=siUNAX=sire=0 0(C+resi)ysiC+resi

  3. Pon . Esta es una buena aproximación inicial con la cual comenzar la iteración para resolver ( A + D ) xx0=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=bBC

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

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 siCC=UNAsisiXXX=siyyX(C+resi)y=si. Por lo tanto, en este caso, el proceso anterior da como punto de partida , que es la solución exacta.X0 0=siy=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.XsiX

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

Arnold Neumaier
fuente
Gracias, profesor Neumaier. Probaré esto. ¿Podrías darme una breve explicación de cómo funciona esto?
Costis
Además, ¿qué pasa si quiero resolver el mismo sistema para muchos vectores RHS diferentes? es decir, , A x 1 = b 1 , A x 2 = b 2 , etc. ¿Hay alguna información que pueda usar de las soluciones anteriores para acelerar las siguientes? UNAX0 0=si0 0UNAX1=si1UNAX2=si2
Costis
@Costis: Una resolución con la misma matriz es solo el caso especial del problema general. Para su primera pregunta, vea la edición. re=0 0
Arnold Neumaier
@Costis: agregué un poco más de detalle al paso 2. - Si escribe la solicitud, envíe una preimpresión de ma.
Arnold Neumaier
(C+resi)ysi