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

9

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 0Ax0=b0A1Ax0

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 1A 0 x 0 M - 1(A+D)x1=b0DA1MMA1A0x0M1 e um método iterativo convergiria em algumas / poucas iterações?

Costis
fonte
Você está começando com um bom pré-condicionador para e deseja saber como atualizá-lo? Qual classificação é a atualização? (A Rank atualização é "pequeno" em comparação com uma matriz de tamanho , mas não pequeno em termos de contagem de iteração.)1000 10 9A1000109
Jed Brown
10 6 10 7A tem o tamanho a e a atualização é <1000 (provavelmente <100) elementos. Estou usando um tipo diagonal de pré-condicionador para A que funciona muito bem, portanto, atualizá-lo seria trivial, mas fiquei pensando se há algo melhor que eu possa fazer em vez de resolver o novo sistema do zero. 106107
Costis
2
A solução de um sistema não diz muito sobre isso. Se você resolver o mesmo sistema várias vezes, o mapa inverso nesses vetores (e / ou espaços de Krylov associados) fornecerá algumas informações que podem ser usadas para acelerar a convergência. Quantos sistemas você está resolvendo em cada caso?
precisa
Atualmente eu só estou resolvendo para um RHS ( vector) com cada matriz antes de modificar . A AbAA
Costis

Respostas:

4
  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 jBCbjcj=Abj

  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 B(A+D)x=bAx=bD=0(C+DB)ybC+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 Cx0=By(A+D)x=bBC

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

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.BCC=ABBxxx=Byyx(C+DB)y=bx0=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.xBx

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

Arnold Neumaier
fonte
Obrigado, Prof. Neumaier. Eu vou tentar isso. Você poderia talvez me dar uma breve explicação sobre como isso funciona?
Costis
Além disso, e se eu quiser resolver o mesmo sistema para muitos vetores RHS diferentes? ie , , , etc. Existe Alguma informação que eu possa usar das soluções anteriores para acelerar as subseqüentes? A x 1 = b 1 A x 2 = b 2Ax0=b0Ax1=b1Ax2=b2
Costis
@ Costis: Uma resolução com a mesma matriz é apenas o caso especial do problema geral. Para sua primeira pergunta, veja a edição. D=0
Arnold Neumaier
@ Costis: adicionei um pouco mais de detalhes à etapa 2. - Se você escrever o aplicativo, envie uma pré-impressão.
Arnold Neumaier
Obrigada pelo esclarecimento! Por que não consigo resolver o sistema sobredeterminado usando uma abordagem baseada em fatoração QR e usando todas as linhas em vez de apenas um subconjunto? Acho que, à medida que o número de colunas de C e B aumenta, talvez eu precise me livrar de algumas linhas para acelerar a operação. Claro, vou escrever uma descrição do sistema e enviar por e-mail para você. Na verdade, acho que é possível criar um esquema específico de aplicativo que funcione melhor do que o caso mais geral. Obrigado! (C+DB)yb
Costis