Mínimos quadrados, de três maneiras

O problema linear de mínimos quadrados procura os coeficientes $\hat\beta$ que minimizam $\lVert X\beta - y \rVert_2$ para uma matriz alta $X \in \mathbb{R}^{m \times p}$. Anulando o gradiente obtêm-se as equações normais,

e em NumPy ocupam uma linha:

beta = np.linalg.solve(X.T @ X, X.T @ y)

Essa linha é curta, rápida e a menos exata dos três métodos habituais. Esta nota mede quanta exatidão se perde, explica porquê e diz quando a perda não importa.

O ajuste polinomial é a maneira clássica de tornar difícil um problema de mínimos quadrados. Em $m = 100$ pontos igualmente espaçados $t_i \in [0, 1]$, as colunas $1, t, t^2, \dots, t^{p-1}$ da base monomial tornam-se cada vez mais parecidas à medida que $p$ cresce, e a matriz aproxima-se da deficiência de característica.

Para isolar o método de resolução, escolhi o segundo membro de forma que a solução exata seja conhecida: $\beta = (1, \dots, 1)$ e $y = X\beta$, com resíduo nulo. Qualquer erro em $\hat\beta$ vem então da aritmética de vírgula flutuante, e não dos dados.

import numpy as np

def errors(p, m=100):
    t = np.linspace(0.0, 1.0, m)
    X = np.vander(t, p, increasing=True)
    beta = np.ones(p)
    y = X @ beta
    rel = lambda b: np.linalg.norm(b - beta) / np.linalg.norm(beta)

    normal = np.linalg.solve(X.T @ X, X.T @ y)
    Q, R = np.linalg.qr(X)
    qr = np.linalg.solve(R, Q.T @ y)
    svd = np.linalg.lstsq(X, y, rcond=None)[0]
    return np.linalg.cond(X), rel(normal), rel(qr), rel(svd)

np.linalg.qr calcula uma fatorização de Householder, e np.linalg.lstsq usa o método do LAPACK baseado na SVD. Os erros relativos nos coeficientes:

Com doze colunas, as equações normais devolvem coeficientes errados no primeiro algarismo, enquanto a QR e a SVD ainda concordam com a solução exata em oito ou nove algarismos. Representando todos os $p$ de 2 a 14 em função do número de condição, surgem duas retas com declives diferentes.

O número de condição $\kappa_2(A) = \sigma_{\max}(A) / \sigma_{\min}(A)$ mede quanto uma matriz pode ampliar perturbações relativas. Formar $X^\top X$ eleva-o ao quadrado.

Um método regressivamente estável para um sistema quadrado com matriz $A$ devolve uma solução com erro relativo da ordem de $\kappa(A)\,u$, em que $u = 2^{-53} \approx 1{,}1 \times 10^{-16}$ é a unidade de arredondamento da precisão dupla. Resolver as equações normais é uma resolução desse tipo com $A = X^\top X$. A QR de Householder e a SVD nunca formam esse produto: trabalham diretamente com $X$ e, num problema de resíduo pequeno, o seu erro é da ordem de $\kappa(X)\,u$.¹

As estimativas são majorantes a menos de constantes modestas, e os erros medidos ficam sobre as guias tracejadas da figura ou abaixo delas. Cada fator de dez em $\kappa(X)$ custa à QR cerca de um algarismo decimal e às equações normais cerca de dois.

Quando o resíduo é grande, até o melhor algoritmo herda um termo proporcional a $\kappa(X)^2$ vezes o tamanho relativo do resíduo. Esse termo pertence ao problema, e não ao método, pelo que nenhuma escolha de algoritmo o elimina.²

Nada disto torna as equações normais erradas. São a ferramenta certa mais vezes do que a tabela sugere.

• A matriz está bem condicionada. Com $\kappa(X) = 100$, as equações normais perdem cerca de quatro algarismos e ainda entregam doze. Preditores padronizados e pouco colineares estão muitas vezes nesta situação.
• Os dados não cabem em memória. $X^\top X$ e $X^\top y$ acumulam-se numa só passagem pelas linhas, em espaço $p \times p$, e resolvem-se no fim.
• A velocidade importa mais do que os últimos algarismos. Para $m \gg p$, formar $X^\top X$ custa cerca de $m p^2$ operações; a QR de Householder custa cerca de $2 m p^2$.

A correção mais barata está muitas vezes no modelo e não no algoritmo. Os monómios $t^k$ em $[0, 1]$ são uma base notoriamente má. Levar os mesmos pontos para $[-1, 1]$ faz descer o número de condição da matriz de doze colunas de $1{,}2 \times 10^{8}$ para $6{,}9 \times 10^{3}$, e uma base de Chebyshev em $[-1, 1]$ leva-o a $2{,}9$. Nesse ponto, os três métodos concordam até ao último algarismo, e as equações normais são tão boas como qualquer outro.

A lição geral é que a exatidão se decide duas vezes: uma quando o problema é formulado, pelo número de condição, e outra quando é resolvido, consoante o algoritmo o eleve ou não ao quadrado.