14. Resolução de sistemas lineares por Quadrados Mínimos

O Problema de Quadrados Mínimos¶

Suponha que na modelagem de um determinado problema chegue-se ao seguinte sistema linear, o qual deseja-se determinar as incógnitas $x$ e $y$:

Esse sistema pode ser reescrito na forma matricial como

Claramente, o sistema não possui solução exata única (a matriz não é sequer quadrada, menos ainda invertível). Analisando do ponto de vista de transformações lineares, sabemos que o sistema será possível e indeterminado se $(1,0,0)\in [(1,3,5),(2,4,6)]$. No entanto,

Logo, os 3 vetores são linearmente independentes e o sistema é impossível, pois não existem vetores $v$ em $\mathbb{R}^{2}$ tais que $Av=(1,0,0)$, onde $A$ é a matriz associada ao sistema.

Supondo que a modelagem foi feita corretamente (sistemas como esses são perfeitamente possíveis e comuns em problemas reais, principalmente quando se tem mais dados observados que incógnitas a serem determinadas), uma solução aproximada seria determinar o vetor pertencente ao espaço gerado pelas colunas de $A$ (nesse caso, uma combinação linear entre os vetores $(1,3,5)$ e $(2,4,6)$) que está “mais próximo” do vetor $(1,0,0)$. Isto é, utilizando a norma vetorial, queremos determinar $\bar{x}\in \mathbb{R}^{2}$ que minimize $\lVert b - A\bar{x} \rVert$ ($b=(1,0,0)$).

Considere, então, o caso geral $Ax=b$ (com $A\in \mathbb{R}^{m\times n}$). O que procuramos é:

Pensando geometricamente, ${} \lVert b-Ax \rVert {}$ corresponde à distância do ponto $b\in \mathbb{R}^{m}$ até o subespaço de $\mathbb{R}^{m}$ gerado pelas colunas de $A$, logo, a menor distância será aquela ortogonal ao subespaço. Isso significa que o vetor $\bar{x}\in \mathbb{R}^{n}$ procurado é tal que o vetor $(b-A\bar{x})\in \mathbb{R}^{m}$ é ortogonal ao espaço coluna de $A$ e, em particular, é ortogonal a cada vetor coluna de $A$.

Portanto, sejam $a_{1},\dots,a_{n}$ as respectivas colunas de $A$, temos:

Utilizando a notação matricial, isso é equivalente a:

O que por sua vez nos dá $A^{T}(b-A\bar{x})=0$. Essa equação, comumente escrita como $\boxed{A^{T}A\bar{x}=A^{T}b}$, é chamada de equação normal e é a base para a solução do problema $Ax=b$ através da abordagem desenvolvida até aqui, denominada de Método dos Quadrados Mínimos (o nome vem da análise do problema do ponto de vista do Cálculo, que procura minimizar o erro quadrático $\lVert b-A\bar{x} \rVert^{2}$, resultando na mesma equação normal).

Soluções da Equação Normal¶

Consideremos primeiramente o caso mais simples, quando $A\in \mathbb{R}^{m\times n}$, com $m\geq n$, possui posto completo (ou seja, $n$). Pelo Teorema 9.12 e a equivalência entre transformações lineares e matrizes, sabemos que o posto de $A^{T}A$, que é $n \times n$, também será $n$. Logo, $A^{T}A$ é invertível e podemos resolver $A^{T}A\bar{x}=A^{T}b$ por:

Apesar de podermos determinar $\bar{x}$ exatamente nesse caso, é importante destacar que ele corresponde à solução da equação normal $A^{T}A\bar{x}=A^{T}b$, e não ao problema original $Ax=b$ que, como vimos, pode nem mesmo possuir soluções. Esse é o caso do sistema linear (2), uma vez que tal $A$ possui 3 linhas e 2 colunas ($m > n$) e posto completo ($\min\{ 3,2 \}=2$), mas o sistema é impossível.

Quando, no entanto, $A$ é quadrada e invertível, a solução da equação normal coincide exatamente com a solução do sistema original, pois:

O mesmo acontece quando $Ax=b$ é possível e indeterminado. Nesse caso, $A^{T}A$ não será invertível e as infinitas soluções de $A^{T}A\bar{x}=A^{T}b$ irão coincidir com as de $Ax=b$.

Uma vez que $Ax=b$ é um sistema impossível, as soluções (ou a solução única, no caso $m \geq n$ e posto completo) de $A^{T}A\bar{x}=A^{T}b$ são soluções aproximadas de $Ax=b$, com base no critério de minimizar $\lVert b - A\bar{x} \rVert$.

Caso geral e pseudoinversa¶

Utilizando a pseudoinversa de Moore-Penrose, dada por Definição 11.9 (Pseudoinversa de Moore-Penrose), podemos resolver $A^{T}Ax=A^{T}b$ no caso geral, que independe da quantidade de linhas e colunas, assim como do posto de $A$.

Concluímos que, quando $m\geq n$ e $A$ tem posto completo, $(A^{T}A)^{-1}A^{T}=A^{+}$. Em geral, calcular a decomposição SVD de $A$ e determinar sua pseudoinversa tem menor custo computacional que calcular $(A^{T}A)^{-1}A^{T}$, uma vez que inversão de matrizes tem alto custo.

Agora, estendemos esse resultado para o caso geral, onde $A^{+}b$ nos dará uma solução única para a equação normal em função de um critério específico.

Com isso, temos uma solução geral para o problema de Quadrados Mínimos, e ainda, única pelo critério de norma mínima (mesmo que a equação normal possua infinitas soluções).

Exemplos numéricos¶

Voltando ao sistema inicial (2), encontraremos sua solução aproximada via quadrados mínimos utilizando a inversa da matriz $A^{T}A$ (pois como já vimos, a matriz associada a esse sistema tem posto completo).

A equação normal associada é dada por:

Então, calculamos a inversa da matriz $A^{T}A$:

Obtendo a solução de quadrados mínimos:

Como constatado no Teorema 14.1, essa solução é a mesma fornecida pela pseudoinversa.

Para ilustrar a solução de quadrados mínimos via pseudoinversa, considere o seguinte sistema, que não possui solução exata e nem posto completo:

Calculando a pseudoinversa, obtemos:

Logo, a solução de quadrados mínimos de menor norma é dada por:

Assim como no tópico de Compressão de Imagens, podemos utilizar a linguagem Python e as funções de Álgebra Linear da biblioteca NumPy, dessa vez para resolver quadrados mínimos via pseudoinversa numericamente. Isso normalmente é feito em problemas reais, cujas matrizes possuem dimensões maiores e o cálculo exato da solução é pouco viável:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
import numpy as np

# Definir as matrizes
A = np.array([[1, 2],
              [2, 4],
              [3, 6]], dtype=float)
              
b = np.array([[1],
              [0],
              [0]], dtype=float)

# Calcular a pseudoinversa 
A_pinv = np.linalg.pinv(A)

# Calcular a solução
x = A_pinv @ b

print("\nPseudoinversa de A:")
print(A_pinv)
print("Solução:")
print(x)

Note que o resultado é determinado numericamente e normalmente não corresponde à solução exata, mas é suficientemente preciso na maioria dos casos práticos.

Obtemos:

Com uma precisão de sete casas decimais em relação à solução exata encontrada.