Resolvendo um sistema com uma pequena atualização diagonal de classificação

Suponha que eu tenha o sistema linear grande e esparso original: . Agora, eu não tenho pois A é muito grande para fatorar ou qualquer tipo de decomposição de , mas suponha que eu tenha a solução encontrada com uma solução iterativa.A - 1 A x$A\textbf{x}_0=\textbf{b}_0$$A^{-1}$$A$$\textbf{x}_0$

Agora, desejo aplicar uma pequena atualização de classificação à diagonal de A (altere algumas das entradas diagonais): que é uma matriz diagonal com principalmente 0 na diagonal e alguns valores diferentes de zero. Se eu tivesse , seria capaz de tirar proveito da fórmula de Woodbury para aplicar uma atualização ao inverso. No entanto, eu não tenho isso disponível. Existe algo que eu possa fazer antes de resolver todo o sistema novamente? Existe alguma maneira, talvez, de conseguir um pré-condicionador que seja fácil \ mais fácil de inverter, como , de modo que tudo o que eu precisaria fazer se tivesse fosse aplicar D A - 1 M M A 1 ≈ A 0 x 0 M - $(A+D)\textbf{x}_1=\textbf{b}_0$$D$$A^{-1}$$M$$MA_1 \approx A_0$$\textbf{x}_0$$M^{-1}$ e um método iterativo convergiria em algumas / poucas iterações?

1. Salve nas colunas das duas matrizes e todos os vetores aos quais você aplicou a matriz nas iterações anteriores e os resultados .C b j c j = A b $B$$C$$b_j$$c_j=Ab_j$

2. Para cada novo sistema (ou , que é o caso especial ), resolva aproximadamente o sistema linear sobredeterminado , por exemplo, selecionando um subconjunto das linhas (possivelmente todas) e usando um método de mínimos quadrados densos. Observe que apenas a parte selecionada do precisa ser montada; então esta é uma operação rápida! A x = b ′ D = 0 ( C + D B ) y ≈ b ′ C + D $(A+D)x'=b'$$Ax=b'$$D=0$$(C+DB)y\approx b'$$C+DB$

3. Coloque . Essa é uma boa aproximação inicial com a qual iniciar a iteração para resolver . Caso outros sistemas precisem ser processados, use os produtos vetoriais de matriz nesta nova iteração para estender as matrizes e no subsistema resultante.( A + D ) x ′ = b ′ B $x_0=By$$(A+D)x'=b'$$B$$C$

Se as matrizes e não couberem na memória principal, armazene no disco e selecione o subconjunto de linhas com antecedência. Isso permite manter no núcleo a parte relevante de e necessária para formar o sistema de mínimos quadrados, e o próximo pode ser calculado por uma passagem por com pouco uso da memória do núcleo.$B$$C$$B$$B$$C$$x_0$$B$

As linhas devem ser selecionadas de forma que correspondam aproximadamente a uma discretização grosseira do problema completo. Tomar cinco vezes mais linhas que o número total de multiplicações de vetores de matriz esperadas deve ser suficiente.

Edit: Por que isso funciona? Por construção, as matrizes e são relacionadas por . Se o subespaço estendido pelas colunas de contiver o vetor exato da solução (uma situação rara, mas simples), então terá a forma durante algum . Substituindo isso na equação que define obtém-se a equação . Portanto, neste caso, o processo acima fornece como ponto de partida , que é a solução exata.$B$$C$$C=AB$$B$$x'$$x'$$x'=By$$y$$x'$$(C+DB)y= b'$$x_0=By=x'$

Em geral, não se pode esperar que esteja no espaço da coluna de , mas o ponto inicial gerado será o ponto nesse espaço de espaço mais próximo de , em uma métrica determinada pelas linhas selecionadas. Portanto, é provável que seja uma aproximação sensata. À medida que mais sistemas são processados, o espaço da coluna aumenta e é provável que a aproximação melhore muito, de modo que se possa esperar convergir em menos e menos iterações.$x'$$B$$x'$

Edit2: Sobre o subespaço gerado: Se alguém resolver cada sistema com um método Krylov, os vetores usados ​​para obter o ponto de partida para o segundo sistema abrangem o subespaço Krylov do primeiro lado direito. Assim, obtém-se uma boa aproximação sempre que esse subespaço de Krylov contiver um vetor próximo à solução do seu segundo sistema. Em geral, os vetores usados ​​para obter o ponto de partida para o sistema st abrangem um espaço que contém o subespaço de Krylov dos primeiros lados do lado direito.$(k+1)$$k$