A solução de um sistema linear de equações pode ser aproximada apenas para as primeiras variáveis?

15

Eu tenho um sistema linear de equações de tamanho mxm, onde m é grande. No entanto, as variáveis ​​nas quais estou interessado são apenas as primeiras n variáveis ​​(n é pequeno comparado a m). Existe uma maneira de aproximar a solução dos primeiros m valores sem precisar resolver todo o sistema? Nesse caso, essa aproximação seria mais rápida do que resolver o sistema linear completo?

Paulo
fonte
2
A menos que sua função forçadora também esteja restrita às primeiras n variáveis. Se for, você pode formar o complemento de Schur, embora seja provavelmente denso. Se o seu operador original for escasso, pode não valer a pena.
22612 Jack Poulson
11
Suponho que você possa usar a eliminação gaussiana a partir do canto inferior direito da matriz. Isso seria ~ 2x mais rápido que a eliminação gaussiana normal se você se importar apenas com os primeiros elementos e parar no meio. Não sei como isso se compara aos métodos iterativos.
Dan
4
@OscarB: Por favor, não. A regra de Cramer é uma atrocidade na aritmética de ponto flutuante. Eu nunca ouvi falar dele sendo usado para cálculos sérios, e é preciso uma quantidade razoável de pensamento para evitar a complexidade fatorial , onde ainda não é competitivo com a eliminação gaussiana.
21412 Jack Poulson
11
@Paul: A maioria das reduções de pedidos de modelos é usada no contexto de grandes sistemas ODE ou DAE. Às vezes, as metodologias de redução são motivadas pelos sistemas ODE ou DAE que surgem da discretização dos PDEs. Não vi redução de modelo usada em equações puramente algébricas. (Se você tiver, envie-me referências, porque estou fazendo minha tese sobre métodos de redução de modelo e ficaria muito interessado em vê-la.) equações algébricas como um caso degenerado de um sistema de equações algébricas diferenciais.
precisa
11
@JackPoulson - você se importa de resumir seu comentário como resposta? Eu acho que é a solução mais correta e não quero que ela se perca nos comentários.
Aron Ahmadia

Respostas:

13

Como outros já apontaram, isso é difícil de resolver com um solucionador direto. Dito isto, não é tão difícil fazer isso com solucionadores iterativos. Para esse fim, observe que a maioria dos solucionadores iterativos, de uma maneira ou de outra, minimiza o erro com relação a alguma norma. Muitas vezes, essa norma é induzida pela própria matriz, mas às vezes também é apenas a norma do vetor l2. Mas isso não precisa ser o caso: você pode escolher em qual norma deseja minimizar o erro (ou residual) e, por exemplo, escolher uma norma na qual pesa os componentes de que gosta 1 e todos os outros com 1e-12, ou seja, por exemplo, algo como (1e-24)N i = 6 x 2 i e produto escalar correspondente. Em seguida, escreva todas as etapas do solucionador iterativo com relação a essa norma e produto escalar, e você obtém um solucionador iterativo que presta muito mais atenção aos elementos vetoriais com os quais você se importa do que aos outros.||x||2=i=15xi2+i=6Nxi2

A questão, é claro, é se você precisa de menos iterações do que com o produto normal / escalar que pesa todos os componentes igualmente. Mas esse deve realmente ser o caso: digamos que você se preocupe apenas com os cinco primeiros elementos do vetor. Em seguida, você deve precisar de no máximo cinco iterações para reduzir o erro em um fator de 1e12, pois cinco iterações são necessárias para o sistema 5x5 que as descreve. Isso não é uma prova, mas tenho certeza de que você deve se livrar de um número muito menor de iterações se o peso na norma (1e-12 acima) for menor que a tolerância com a qual você deseja resolver o sistema linear iterativamente .

Wolfgang Bangerth
fonte
2
Hmm, bom ponto. Eu estaria interessado em ver um exemplo real, pois me preocupo um pouco com os efeitos de apenas tentar resolver alguns graus de liberdade; mesmo que o resíduo possa ser pequeno, talvez a norma do erro ainda seja bastante grande (faça isso para ignorar efetivamente a maioria do operador).
21412 Jack Poulson
Intuitivamente, isso só parece funcionar se os componentes do sistema muito pequeno realmente dominam a resposta em um sentido L2 (ou a norma em que você entende que seu erro é medido). Caso contrário, acho que a preocupação de Jack é válida, mas eu definitivamente estaria interessado em ver uma prova numérica disso ...
Aron Ahmadia
Seria necessário ter um método que minimizasse o erro , não o residual. Eu acho que o MinErr pode ser um bom ponto de partida.
Wolfgang Bangerth
@WolfgangBangerth: Não estou familiarizado com o MINERR: esta é a principal referência?
21412 Jack Poulson
11
Mesmo isso não basta, porque você será impreciso. Você não pode obter alguns componentes com precisão usando essa ponderação.
precisa saber é o seguinte
17

Formando o complemento Schur

Suponha que você tenha permutado e particionado sua matriz no formato

A=(A11A12A21A22),

tal que A22 contenha seus graus de liberdade de interesse e seja muito menor que , é possível formar o complemento de SchurA11

S22:=A22A21A111A12,

através de uma fatoração LU de aparência parcial ou da fórmula explícita e, em seguida, pode ser entendido no seguinte sentido:S22

S22x=y(A11A12A21A22)(x)=(0y),

onde representa a parte 'desinteressante' da solução. Assim, desde que o lado direito seja diferente de zero nos graus de liberdade do complemento Schur S 22 , precisamos apenas resolver contra S 22 para obter a parte da solução correspondente a esses graus de liberdade.S22S22

Complexidade computacional em caso denso não estruturado

Configuração para a altura de um e n para a altura de Um 22 , então o método padrão para calcular S 22 é primeiro factor G 11 L 11 : = A 11 (Vamos ignorar pivotante por agora) em aproximadamente 2 / 3 ( N - nNAnA22S22L11U11:=A11 trabalhos, depois formar2/3(Nn)3

S22:=A22(A21U111)(L111A12)=A22A21A111A12

usando duas soluções triangulares exigindo que trabalhem cada uma e, em seguida, executando a atualização para A 22 em 2 n 2 ( Nn(Nn)2A22 .2n2(Nn)

Assim, o trabalho total é de cerca de . Quando n é muito pequena, N - N N , de modo que o custo pode ser visto como cerca de 2 / 3 N 3 , que é o custo de um completo fatoração.2/3(Nn)3+2n(Nn)2+2n2(Nn)nNnN2/3N3

O benefício é que, se houver um número muito grande de lados do lado direito a ser resolvido com o mesmo sistema de equações, o poderá ser reutilizado potencialmente várias vezes, onde cada solução exigiria apenas 2 n 2 de trabalho (em vez de 2 N 2 funcionar) se S 22 for fatorado.S222n22N2S22

Complexidade computacional no caso escasso (típico)

Se seu sistema esparso surgir de algum tipo de aproximação de diferença finita ou elemento finito, os solucionadores diretos esparsos quase certamente serão capazes de explorar parte da estrutura; Sistemas 2d pode ser resolvido com de trabalho e O ( N log N ) de armazenamento, enquanto que os sistemas 3D pode ser resolvido com S ( N 2 ) trabalho e S ( N 4 /O(N3/2)O(NlogN)O(N2)de armazenamento. Os sistemas fatorados podem ser resolvidos com a mesma quantidade de trabalho que os requisitos de armazenamento.O(N4/3)

O objetivo de trazer à tona as complexidades computacionais é que, se e você tem um sistema 2d, então, como o complemento Schur provavelmente será denso, a complexidade da solução dada ao complemento fatorado Schur seráO(n2)=O(N), que está faltando apenas um fator logarítmico em vez de resolver o problema completo sistema! Em 3d, que exigeO(N)de trabalho em vez deS(N 4 / 3nNO(n2)=O(N)O(N) .O(N4/3)

Portanto, é importante ter em mente que, no seu caso, onde , haverá apenas economias significativas se você estiver trabalhando em várias dimensões e tiver muitos lados do lado direito para resolver.n=N

Jack Poulson
fonte
11
Este é um ótimo resumo do método do complemento schur e quando é computacionalmente eficiente usá-lo!
Paul
6

A abordagem de redução de modelo

Como Paul perguntou, vou falar sobre o que acontece se você usar métodos de redução de modelo baseados em projeção sobre esse problema. Suponha que você possa criar um projetor modo que o intervalo de P , denotado R ( P ) , contenha a solução para seu sistema linear A x = b e tenha a dimensão k , em que k é o número de incógnitas pelas quais você deseja resolver em um sistema linear.PPR(P)Ax=bkk

Uma decomposição de valor singular de produzirá a seguinte matriz particionada:P

P=[V][diag(1k)000][WT].

As matrizes obscurecidas pelas estrelas são importantes para outras coisas (como estimar erros, etc.), mas, por enquanto, evitaremos lidar com detalhes estranhos. Segue que

P=VWT

um ponto de decomposição completa de .P

Essencialmente, você resolverá o sistema

PAx=Pb

de uma forma inteligente, porque e W também tem a propriedade de que W t V = I . Multiplicando ambos os lados de P Um x = P b por W T e deixando y = V x ser uma aproximação para x rendimentosVWWTV=IPAx=PbWTy=Vx^x

WTAx^=WTb.

Resolva para x , premultiply-lo por V , e você tem y , a sua aproximação de x .x^Vyx

Por que a abordagem do complemento Schur é provavelmente melhor

Para começar, você deve escolher alguma forma. Se a solução para A x = b estiver em R ( P ) , então y = x , y não é uma aproximação. Caso contrário, yx , e você introduzir algum erro de aproximação. Essa abordagem realmente não aproveita toda a estrutura que você mencionou que deseja explorar. Se escolhermos P de modo que seu intervalo seja a base da unidade padrão nas coordenadas de x que você deseja calcular, as coordenadas correspondentes de y terão erros nelas. Não está claro como você gostaria de escolherPAx=bR(P)y=xyyxPxy . Você pode usar um SVD de A , por exemplo, e selecionar P para ser o produto dos primeiros k vetores singulares à esquerda de A e ao lado dos primeiros k vetores singulares à direita de A , assumindo que vetores singulares sejam organizados em ordem decrescente de valor singular. Essa escolha do projetor seria equivalente à realização de decomposição ortogonal adequada em A e minimizaria oerroL 2 na solução aproximada.PAPkAkAA2

Além disso a introdução de erros de aproximação, esta abordagem também apresenta três multiplicações de matriz extra no topo do resolver linear do sistema menor e o trabalho necessários para calcular , e W . A menos que você esteja resolvendo muito o mesmo sistema linear, mudando apenas o lado direito e P ainda seja uma matriz de projeção "boa" para todos esses sistemas, esses custos extras provavelmente tornarão mais caro a solução do sistema reduzido do que a solução do seu sistema. sistema original.VWP

As desvantagens são muito parecidas com a abordagem de JackPoulson, exceto que você não está aproveitando bastante a estrutura que mencionou.

Geoff Oxberry
fonte
4

A resposta longa é ... mais ou menos.

Você pode reorganizar seu sistema de equações de modo que as colunas mais à direita sejam as variáveis ​​que você deseja resolver.k

Etapa 1: Execute a Eliminação Gaussiana para que a matriz seja triangular superior. Etapa 2: resolva por substituição inversa apenas as primeiras (últimas) variáveis ​​nas quais você está interessadok

Isso poupará a complexidade computacional de ter que resolver as últimas variáveis ​​via substituição traseira, o que poderia valer a pena se n for tão grande quanto você diz. Lembre-se de que ainda será necessário realizar uma boa quantidade de trabalho para a etapa 1.nkn

Além disso, mantenha em mente que a restrição da ordem em que você está indo para executar back-substituion pode restringir a forma da matriz (que tira a capacidade de colunas de troca) que poderia possivelmente levar a um sistema mal-condicionado, mas eu não sou tenho certeza disso - apenas algo a ter em mente.

drjrm3
fonte
O(n3)O(n2)n
é por isso que a resposta é "espécie de" em vez de "sim" =)
drjrm3
Faz sentido que isso possa ser feito dessa maneira ... No entanto, a maior parte do cálculo em uma Eliminação Gaussiana está na fase de eliminação direta, produzindo uma complexidade de O (n ^ 3) apesar da fase de substituição reversa truncada. Eu estava esperando que houvesse um método mais rápido ...
Paul