Quando as transformações ortogonais superam a eliminação gaussiana?

Como sabemos, os métodos de transformações ortogonais (rotações de Givens e reflexões de Housholder) para sistemas de equações lineares são mais caros que a eliminação gaussiana, mas teoricamente possuem propriedades de estabilidade mais agradáveis, no sentido de que não alteram o número de condição do sistema. Embora eu conheça apenas um exemplo acadêmico de uma matriz que é estragada pela eliminação gaussiana com rotação parcial. E existe uma opinião comum de que é muito improvável encontrar esse tipo de comportamento na prática (consulte as notas desta palestra [pdf] ).

Então, onde procuraremos a resposta sobre o assunto? Implementações paralelas? Atualizando? ..

Precisão

Trefethen e Schreiber escreveram um excelente artigo, Estabilidade média da eliminação gaussiana , que discute o lado preciso da sua pergunta. Aqui estão algumas de suas conclusões:

1. "Para a fatoração QR com ou sem pivô da coluna, o elemento máximo médio da matriz residual é , enquanto na eliminação gaussiana é . Essa comparação revela que a eliminação gaussiana é levemente instável , mas a instabilidade só seria detectável para problemas de matriz muito grandes resolvidos com baixa precisão. Para a maioria dos problemas práticos, a eliminação gaussiana é altamente estável, em média. "(Ênfase minha)S ( n $O(n^{1/2})$$O(n)$

2. "Após os primeiros passos da eliminação gaussiana, os elementos remanescentes da matriz são aproximadamente normalmente distribuídos, independentemente de terem começado dessa maneira".

Há muito mais no artigo que não posso capturar aqui, incluindo a discussão da matriz de pior caso que você mencionou, por isso recomendo fortemente que você a leia.

atuação

Para matrizes reais quadradas, a LU com pivotagem parcial requer aproximadamente flops, enquanto o QR baseado em chefe de família requer aproximadamente flops. Assim, para matrizes quadradas razoavelmente grandes, a fatoração QR será apenas duas vezes mais cara que a fatoração LU. 4 / 3 n $2/3 n^3$$4/3 n^3$

Para matrizes, em que , LU com pivô parcial requer 3/3 flops, versus do QR (que ainda é o dobro da fatoração da LU) . No entanto , é surpreendentemente comum que aplicações produzam matrizes magras muito altas ( ), e Demmel et al. tenha um bom artigo, fatoração QR paralela e seqüencial para evitar a comunicação , que (na seção 4) discute um algoritmo inteligente que exige apenas que mensagens sejam enviadas quando os processadores são usados, versus as mensagens das abordagens tradicionais . A despesa é quem ≥ n m n 2 - N 3 / 3 2 m n 2 - 2 n 3 / 3 m » n log p p n log P ó ( n 3 log P ) $m \times n$$m \ge n$$mn^2 - n^3/3$$2mn^2 - 2n^3/3$$m \gg n$$\log p$$p$$n \log p$$O(n^3 \log p)$ flops extras são executados, mas para muito pequeno isso costuma ser preferido ao custo de latência do envio de mais mensagens (pelo menos quando apenas uma única fatoração QR precisa ser executada).$n$

Estou surpreso que ninguém tenha mencionado problemas lineares mínimos quadrados , que ocorrem com freqüência na computação científica. Se você deseja usar a eliminação gaussiana, deve formar e resolver as equações normais, que se parecem com:

onde é uma matriz de pontos de dados que corresponde a observações de variáveis ​​independentes, é um vetor de parâmetros a serem encontrados é um vetor de pontos de dados que corresponde a observações de uma variável dependente.x $A$$x$$b$

Como Jack Poulson frequentemente indica, o número de condição de é o quadrado do número de condição de , de modo que as equações normais podem ser desastrosamente mal condicionadas. Nesses casos, embora as abordagens baseadas em QR e SVD sejam mais lentas, elas produzem resultados muito mais precisos.$A^{T}A$$A$

Como você mede o desempenho? Rapidez? Precisão? Estabilidade? Um teste rápido no Matlab fornece o seguinte:

>> N = 100;
>> A = randn(N); b = randn(N,1);
>> tic, for k=1:10000, [L,U,p] = lu(A,'vector'); x = U\(L\b(p)); end; norm(A*x-b), toc
ans =
   1.4303e-13
Elapsed time is 2.232487 seconds.
>> tic, for k=1:10000, [Q,R] = qr(A); x = R\(Q'*b); end; norm(A*x-b), toc             
ans =
   5.0311e-14
Elapsed time is 7.563242 seconds.

Portanto, resolver um único sistema com uma decomposição de LU é cerca de três vezes mais rápido que resolvê-lo com uma decomposição de QR, ao custo de meio dígito decimal de precisão (neste exemplo!).

O artigo que você cita defende a Eliminação Gaussiana dizendo que, embora seja numericamente instável, tende a se dar bem em matrizes aleatórias e, como a maioria das matrizes que se pode pensar são como matrizes aleatórias, devemos ficar bem. Essa mesma afirmação pode ser dita sobre muitos métodos numericamente instáveis.

Considere o espaço de todas as matrizes. Esses métodos funcionam bem em quase todos os lugares. Isso é 99,999 ...% de todas as matrizes que se poderia criar não terá problemas com métodos instáveis. Existe apenas uma fração muito pequena de matrizes para a qual a GE e outras terão dificuldades.

Os problemas com que os pesquisadores se preocupam tendem a estar nessa pequena fração.

Não construímos matrizes aleatoriamente. Construímos matrizes com propriedades muito especiais que correspondem a sistemas não-aleatórios muito especiais. Essas matrizes geralmente são mal condicionadas.

Geometricamente, você pode considerar o espaço linear de todas as matrizes. Há um subespaço de volume / medida zero de matrizes singulares que corta esse espaço. Muitos problemas que construímos estão agrupados em torno desse subespaço. Eles não são distribuídos aleatoriamente.

Como exemplo, considere a equação ou dispersão do calor. Esses sistemas tendem a remover informações do sistema (todos os estados iniciais gravitam para um único estado final) e, como resultado, as matrizes que descrevem essas equações são enormemente singulares. Esse processo é muito improvável em uma situação aleatória, mas onipresente em sistemas físicos.