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?
linear-algebra
approximation
Paulo
fonte
fonte
Respostas:
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=∑5i=1x2i+ ∑Ni=6x2i
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 .
fonte
Formando o complemento Schur
Suponha que você tenha permutado e particionado sua matriz no formato
tal queA22 contenha seus graus de liberdade de interesse e seja muito menor que , é possível formar o complemento de SchurA11
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
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.⋆ S22 S22
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 - nN A n A22 S22 L11U11:=A11 trabalhos, depois formar2/3(N−n)3
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(N−n)2 A22 .2n2(N−n)
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(N−n)3+2n(N−n)2+2n2(N−n) n N−n≈N 2/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.S22 2n2 2N2 S22
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 / 3n≈N−−√ O(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−−√
fonte
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.P P R(P) Ax=b k k
Uma decomposição de valor singular de produzirá a seguinte matriz particionada:P
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
um ponto de decomposição completa de .P
Essencialmente, você resolverá o sistema
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 rendimentosV W WTV=I PAx=Pb WT y=Vxˆ x
Resolva para x , premultiply-lo por V , e você tem y , a sua aproximação de x .xˆ V y x
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, y ≠ x , 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 escolherP Ax=b R(P) y=x y y≠x P x y . 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.P A P k A k A A 2
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.V W P
As desvantagens são muito parecidas com a abordagem de JackPoulson, exceto que você não está aproveitando bastante a estrutura que mencionou.
fonte
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.n−k n
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.
fonte