Condicionamento da matriz
Arredondamento do float
Um exemplo em que operações entre floats se tornam um grande problema é o cancelamento catastrófico:
Suponhamos que temos o número 1000,1 onde 0,1 é erro de arredondamento, até aqui temos o erro absoluto (0,1) e o erro relativo (0,01%) pequenos.
Tivemos que subtrair por algum número próximo, por exemplo 1000,1 - 999 = 1,1.
Note que agora nosso erro absoluto (0,1) continua pequeno, porém nosso erro relativo (10%) está grande.
Tivemos que multiplicar esse número, por exemplo 1,1 x 200 = 220.
Note que agora nosso erro absoluto (20) e o erro relativo (10%) estão grandes, basicamente amplificamos o erro.
Um resultado desse exemplo é que a quantidade de operações entre floats não implica instabilidade, quais operações estão sendo feitas podem ter mais influência, poderíamos somar 1,0 um milhão de vezes e os erros estariam estáveis.
Número de Condicionamento
Queremos uma forma de dizer se os elementos da matriz estão suscetíveis aos erros de arredondamento.
Exemplo:
Ao calcularmos a inversa:
Relembre
Note que aparentemente há uma diferença na grandeza dos números, isso não é uma coincidência, no processo de Eliminação de Gauss-Jordan para acharmos a Inversa chegamos a:
Para termos o pivô com valor 1, temos que multiplicar a segunda linha por 10.000. Se temos linhas quase Linearmente Dependentes, então teremos mais chance desse evento ocorrer.
Uma consequência disso é que no sistema \(A\mathbf{x}=\mathbf{b}\), pequenas mudanças em \(\mathbf{b}\) causam grandes alterações na solução \(\mathbf{x}\).
Exemplo
Note que uma pequena alteração de -0.0001 no elemento \(A_{2x2}\) torna a matriz não inversível.
Algumas matrizes são mais suscetíveis a ter esses erros do que outras, então queremos alguma forma de medir essa propriedade da matriz que chamaremos de Condicionamento e a notação será \(\kappa (A)\).
Queremos saber se podemos dizer algo do "erro relativo de \(\mathbf{x}\)", dado que temos o "erro relativo de \(\mathbf{b}\)" e o "Condicionamento de A". Algo como \(\frac{||\mathbf{\Delta x}||}{||\mathbf{x}||}\) e \(\frac{||\mathbf{\Delta b}||}{||\mathbf{b}||} \kappa (A)\).
Um resultado importante é que conseguimos achar uma relação: \(\frac{||\mathbf{\Delta b}||}{||\mathbf{b}||} \kappa (A)\) é o limite superior para o erro relativo de \(\mathbf{x}\), onde \(\kappa (A) = ||A|| \cdot ||A^{-1}||\) (note como o "tamanho" dos elementos da matriz influenciam como no exemplo).
Demonstração
Usaremos a inequação de normas matriciais:
O Sistema Linear na prática prática pode ser visto como:
Note que podemos separar em 2 Sistemas Lineares:
O primeiro sistema pode dar os denominadores da relação que queremos, o segundo os numeradores:
Para chegarmos à relação do erro relativo juntamos as duas inequações:
Relembre
O número de condicionamento pode ser escrito como \(\frac{\sigma_{max}}{\sigma_{min}}\), então o \(\sigma_{min}\) muito pequeno pode causar o número de condicionamento maior (quanto maior, pior a estabilidade).
Como vimos no exemplo, linhas quase Linearmente Dependentes aumentam o número de condicionamento, mas colunas quase Linearmente Dependentes também, isso ocorre porque ambas implicam que existe um Valor Singular muito pequeno. Note que \(A\) e \(A^T\) têm os mesmos Valores Singulares quando formos fazer a Decomposição de Valor Singular, então linhas e colunas podem ser tratadas iguais nesse contexto.
No exemplo dado, note que as colunas são quase paralelas também, se multiplicarmos A por \(\begin{pmatrix}-1 \\ 1\end{pmatrix}\) (usando a Combinação Linear das colunas) o resultado será quase o vetor nulo, se normalizarmos esse vetor entrada, então a norma do resultado pode ser o \(\sigma_{min}\).
A intuição geométrica é que uma matriz mal-condicionada será um elipsoide quase que completamente achatado, então nas direções achatadas a densidade de pontos seria muito grande, teríamos que ter uma precisão computacionalmente inviável para conseguir distinguir os pontos e achar uma aproximação melhor para \(\mathbf{x}\).
Resolução de Sistemas Lineares
Caso o Condicionamento da matriz for baixo, então podemos resolver o sistema achando a inversa e isolando \(\mathbf{x}\). Porém, caso o condicionamento da matriz for alto, as operações da Eliminação de Gauss-Jordan para achar a inversa são de risco, e a multiplicação final \(A^{-1}\mathbf{b}\) amplifica os erros de todos os elementos.
Queremos evitar esse acúmulo de erros, para isso podemos fazer apenas a Eliminação de Gauss, pois tem mais ou menos a metade de operações elementares e evita a operação de tornar o pivô 1.
Lembre que sempre podemos computar a solução de um sistema escalonado (triangular em geral) com uma fórmula pronta sem usar inversa.
Exemplo
Se fizermos a multiplicação do produto interno, sempre teremos uma forma de isolar o \(x_n\), pois poderemos usar os resultados anteriores:
A=LU
Poderíamos ter feito o escalonamento na matriz aumentada, porém se fizermos isso não poderemos reutilizar a computação para um \(\mathbf{b}\) diferente.
Como não alteramos \(\mathbf{b}\), então não podemos apenas substituir A pela sua forma escalonada, devemos usar a decomposição LU.
Relembre
Uma propriedade importante é uma matriz triangular tende a ter o número de condicionamento baixo, note que dificilmente teremos linhas ou colunas quase Linearmente Dependentes, assim evitamos amplificar os erros da Eliminação de Gauss em uma multiplicação com um vetor.
Essa forma de resolver um sistema linear tem mais ou menos 1/4 das operações entre floats da resolução por Inversa, porém ainda estamos na ordem de grandeza cúbica.
Exemplo
Complexidade da Eliminação de Gauss-Jordan em uma matriz nxn:
-
teremos n pivôs nessa matriz
-
para cada pivô teremos que zerar os n-1 elementos da coluna do pivô com uma operação elementar (operação na linha específica)
-
cada operação elementar tem n operações
Então no total teremos \(n(n-1)n=n^3-n^2\) operações, podemos abstrair para \(O(n^3)\).
A=QR
Outra forma de resolver Sistemas Lineares com matrizes triangulares é usando a decomposição QR.
Note como Q é uma matriz de rotação, então seus Valores Singulares são iguais a 1 e \(\kappa (Q)=1\).
Relembre
Para maximizar a estabilidade usamos o método de Householder para a decomposição, pois as matrizes que fazem as reflexões de Householder são ortogonais e têm \(\kappa (Q)=1\).
Relembre
Então essa resolução de Sistema Linear será mais estável que a LU, porém terá mais ou menos 2 vezes mais operações entre floats. Na prática a maioria dos softwares usam LU, pois os piores casos são raros.
Cholesky
Relembre
Outra forma de resolver Sistemas Lineares com matrizes triangulares é usando a decomposição Cholesky, porém só no caso especial que A é Simétrica Positiva Definida.
Como A não precisa de trocas de linhas na Eliminação de Gauss para achar um pivô estável, então temos mais ou menos a metade das operações da resolução por A=LU e evitamos os piores casos de instabilidade.