Álgebra Linear Numérica

Álgebra linear numérica em poucas palavras

A álgebra linear numérica estuda como fazer computações com matrizes e vetores de forma confiável e eficiente em hardware de precisão finita. Na matemática exata, operações como resolver (Ax=b), calcular uma DVS (SVD), ou multiplicar matrizes são bem definidas. Em máquinas reais, você precisa lidar com:

• Erro de arredondamento (a aritmética de ponto flutuante é aproximada)
• Condicionamento (alguns problemas são inerentemente sensíveis a pequenas perturbações)
• Estabilidade (alguns algoritmos amplificam erros muito mais do que outros)
• Escala (modelos e conjuntos de dados grandes tornam métodos “de livro-texto” lentos demais ou exigentes demais em memória)

Esses problemas aparecem o tempo todo em AM (ML): mínimos quadrados, ACP (PCA), branqueamento (whitening), estimação de covariância, otimização de segunda ordem, camadas implícitas, verossimilhanças gaussianas e operações em lote grandes em kernels de GPU.

Este artigo conecta a teoria (condicionamento e estabilidade) à orientação prática: quais decomposições usar, o que evitar (por exemplo, inversas explícitas), como escolhas de precisão afetam resultados e como diagnosticar falhas.

Contexto relacionado: Vetores, Matrizes, Tensores e Autovalores & DVS.

Aritmética de ponto flutuante: onde “números reais” se tornam aproximados

A maior parte do código de AM usa tipos IEEE-754 de ponto flutuante:

• float64 (double): ~16 dígitos decimais de precisão
• float32 (single): ~7 dígitos decimais de precisão (padrão comum em GPUs)
• float16 / bfloat16: ~3–4 dígitos decimais de precisão (comum para treinar em escala)

Um modelo mental útil é:

[ \mathrm{fl}(a \ \text{op} \ b) = (a \ \text{op} \ b)(1 + \delta), \quad |\delta| \lesssim u ]

onde (u) é o épsilon de máquina (machine epsilon) (aproximadamente (10^{-16}) para float64, (10^{-7}) para float32).

Modos comuns de falha numérica

• Cancelamento catastrófico: subtrair números quase iguais perde dígitos significativos.
  • Exemplo: (1.0000001 - 1.0000000) em float32 mantém apenas alguns dígitos relevantes.
• Estouro/subfluxo (overflow/underflow): valores excedem o intervalo representável ou “escorrem” para zero.
• Não associatividade: ((a+b)+c \neq a+(b+c)) devido ao arredondamento; reduções paralelas em GPUs podem alterar resultados.
• Erro de acumulação: somar muitos valores pequenos (por exemplo, produtos escalares) aumenta o erro com a dimensão.

Em hardware moderno, a multiplicação e soma fundidas (FMA fused multiply-add) frequentemente melhora a precisão de produtos escalares, mas a ordem de redução ainda importa.

Condicionamento: sensibilidade do *problema*

Condicionamento responde: Se a entrada é perturbada levemente (por ruído, arredondamento, incerteza nos dados), o quanto a saída pode mudar? Mesmo um algoritmo perfeito não consegue superar um problema fundamentalmente mal condicionado.

Número de condição para resolver sistemas lineares

Para resolver (Ax=b), uma medida clássica é o número de condição:

[ \kappa(A) = |A| \cdot |A^{-1}| ]

Frequentemente computado na norma-2 como:

[ \kappa_2(A) = \frac{\sigma_{\max}(A)}{\sigma_{\min}(A)} ]

onde (\sigma_{\max}, \sigma_{\min}) são valores singulares (ver Autovalores & DVS).

Interpretação:

• (\kappa(A) \approx 1): bem condicionado (estável, confiável)
• (\kappa(A) \gg 1): mal condicionado (pequenas perturbações em (A) ou (b) podem causar grandes mudanças em (x))
• Se (\sigma_{\min}(A)) estiver próximo de 0, (A) é quase singular e a solução é extremamente sensível.

Um limite padrão (informal, mas prático) é:

[ \frac{|\Delta x|}{|x|} \lesssim \kappa(A)\left(\frac{|\Delta A|}{|A|} + \frac{|\Delta b|}{|b|}\right) ]

Então, se (\kappa(A)) é (10^8), mesmo o arredondamento de float32 ((\sim 10^{-7})) pode produzir erro relativo de ordem 1.

Exemplo prático: características mal condicionadas

Em AM, o mau condicionamento muitas vezes vem de características altamente correlacionadas, escalonamento ruim ou parâmetros redundantes.

import numpy as np

np.random.seed(0)
n = 200
x1 = np.random.randn(n)
x2 = x1 + 1e-6*np.random.randn(n)  # nearly collinear feature
A = np.c_[x1, x2]
cond = np.linalg.cond(A)
cond

Em geral, você verá um número de condição grande. Isso pode tornar ajustes de mínimos quadrados instáveis, a menos que você regularize (ridge) ou use resolvedores numericamente robustos.

Estabilidade: sensibilidade do *algoritmo*

Estabilidade pergunta: O algoritmo amplifica erros de arredondamento?

Um conceito central é a estabilidade retroativa (backward stability): um algoritmo é retroativamente estável se o resultado computado é a solução exata de um problema próximo.

Para resolver (Ax=b), um algoritmo retroativamente estável produz (\hat{x}) tal que:

[ (A + \Delta A)\hat{x} = b + \Delta b ]

com (|\Delta A|/|A|) e (|\Delta b|/|b|) da ordem do épsilon de máquina.

Por que isso importa:

• Se seu algoritmo é retroativamente estável, então a acurácia final é controlada principalmente pelo condicionamento (que às vezes dá para melhorar via escalonamento ou regularização).
• Se o algoritmo é instável, você pode obter erros grandes mesmo para problemas bem condicionados.

Erro direto vs erro retroativo

• Erro direto (forward error): quão próximo (\hat{x}) está de (x)
• Erro retroativo (backward error): o quanto você precisaria perturbar as entradas para tornar (\hat{x}) exato

O erro retroativo costuma ser mais fácil de controlar e analisar; o erro direto então decorre via condicionamento.

Tarefas canônicas e o que é numericamente seguro

1) Resolver sistemas lineares: nunca forme a inversa

Uma armadilha recorrente é calcular (x = A^{-1}b) formando explicitamente (A^{-1}). Isso geralmente é:

• mais lento
• menos preciso
• menos estável
• mais intensivo em memória

Use um resolvedor dedicado no lugar.

import numpy as np

A = np.random.randn(100, 100)
b = np.random.randn(100)

x_solve = np.linalg.solve(A, b)        # preferred
x_inv   = np.linalg.inv(A) @ b         # avoid

res_solve = np.linalg.norm(A@x_solve - b)
res_inv   = np.linalg.norm(A@x_inv   - b)

res_solve, res_inv

Em frameworks de aprendizado profundo, prefira torch.linalg.solve ou soluções baseadas em Cholesky para matrizes SPD.

2) Mínimos quadrados: evite as equações normais se você se importa com precisão

Mínimos quadrados resolve (\min_x |Ax-b|_2). Uma derivação tentadora leva às equações normais (normal equations):

[ A^\top A x = A^\top b ]

O porém: (\kappa(A^\top A) \approx \kappa(A)^2). Elevar ao quadrado o número de condição pode destruir a acurácia.

Opções melhores:

• Decomposição QR (frequentemente a abordagem estável padrão)
• DVS (melhor para deficiência de posto e diagnósticos, mais custosa)

import numpy as np

np.random.seed(0)
m, n = 200, 20
A = np.random.randn(m, n)
A[:, 0] = A[:, 1] + 1e-8*np.random.randn(m)  # inject near dependency
b = np.random.randn(m)

# Normal equations (unstable when A is ill-conditioned)
x_ne = np.linalg.solve(A.T @ A, A.T @ b)

# Stable least squares
x_ls, *_ = np.linalg.lstsq(A, b, rcond=None)

print("Residual normal eq:", np.linalg.norm(A@x_ne - b))
print("Residual lstsq    :", np.linalg.norm(A@x_ls - b))
print("cond(A):", np.linalg.cond(A))
print("cond(A^T A):", np.linalg.cond(A.T@A))

Em AM, isso aparece em regressão linear, algumas atualizações em forma fechada e métodos clássicos dentro de pipelines. Se você precisar usar equações normais (por exemplo, explorando estrutura), adicione regularização e considere maior precisão.

3) Sistemas simétricos definidos positivos (SPD): use Cholesky, mas observe a definitude

Para matrizes SPD (comuns: matrizes de covariância, matrizes de kernel com jitter, gaussianas), Cholesky é rápida e estável:

[ A = LL^\top ]

Então resolva via substituição para frente/para trás.

Armadilha: uma matriz que teoricamente é SPD pode ser numericamente indefinida devido a arredondamento ou aproximações de modelagem. Uma correção comum é adicionar jitter / regularização de Tikhonov (Tikhonov regularization):

[ A \leftarrow A + \lambda I ]

Isso é comum em regressão por processos gaussianos e ao computar log-determinantes para verossimilhanças gaussianas.

4) Computações de autovalores/DVS: estáveis, mas a interpretação importa

Autovalores e valores singulares são centrais para ACP, condicionamento e estrutura de baixa ordem. Numericamente:

• DVS é, em geral, muito robusta.
• Decomposição em autovalores para matrizes simétricas também é bem estudada e estável.
• Para problemas de autovalores não simétricos, a sensibilidade pode ser pior (autovalores podem ser altamente não robustos).

Em AM em grande escala, você frequentemente usa métodos aproximados (por exemplo, DVS randomizada) para extrair os principais componentes com eficiência (ver Autovalores & DVS).

Condicionamento encontra otimização: por que “mal condicionado” torna o aprendizado mais lento

Mau condicionamento não é apenas sobre resolver (Ax=b). Ele afeta diretamente a dinâmica de otimização.

Para um objetivo quadrático (f(x) = \frac{1}{2}x^\top H x - b^\top x), métodos baseados em gradiente dependem dos autovalores da Hessiana (H). A razão (\lambda_{\max}/\lambda_{\min}) é um número de condição que controla quão “esticada” é a paisagem.

Consequências práticas:

• Descida do Gradiente converge lentamente quando o problema é mal condicionado.
• Pré-condicionamento (explicitamente ou implicitamente) melhora a convergência.
• Camadas de normalização e otimizadores adaptativos podem ser vistos como tentativas parciais de lidar com escalonamento/condicionamento (embora não substituam métodos numéricos adequados em soluções lineares).

Álgebra linear numérica em larga escala em sistemas de AM

Precisão mista: velocidade vs estabilidade

O treinamento moderno frequentemente usa float16/bfloat16 para throughput, com acumulação em float32 por segurança.

Questões-chave:

• Subfluxo (underflow): gradientes podem virar zeros em float16.
• Estouro (overflow): ativações/gradientes podem virar inf/NaN.
• Escalonamento da perda (loss scaling): multiplica a perda para manter gradientes dentro do intervalo representável (comum em treinamento com float16).

Regra prática:

• Faça reduções sensíveis (somas, produtos escalares), estatísticas de normalização e fatorações em float32 ou float64 quando necessário.
• Tenha cautela ao computar (A^\top A) ou matrizes de covariância em baixa precisão; o efeito de elevação ao quadrado pode ser brutal.

Métodos esparsos e iterativos: quando as matrizes não cabem

Em escala, pode ser impossível fatorar (A) (memória/tempo). Para sistemas esparsos grandes, use resolvedores iterativos:

• Gradientes Conjugados (CG) para SPD
• GMRES, BiCGSTAB para sistemas mais gerais

Métodos iterativos dependem fortemente de pré-condicionamento (preconditioning) (transformar o sistema para reduzir o número de condição efetivo). Em AM, ideias de pré-condicionamento aparecem em aproximações de segunda ordem e métodos do tipo gradiente natural.

Álgebra linear em lote em GPUs: não determinismo e reprodutibilidade

Kernels de GPU frequentemente paralelizam reduções de maneiras que não são determinísticas bit a bit. Mesmo que cada operação siga IEEE, ordens diferentes de soma produzem resultados ligeiramente diferentes.

Implicações:

• Pequenas diferenças numéricas podem mudar trajetórias de treinamento (especialmente em regimes caóticos).
• Reprodutibilidade pode exigir modos determinísticos, a um custo de desempenho.
• Sempre avalie correção usando tolerâncias, não igualdade exata.

Armadilhas práticas e “faça isto em vez disso”

Armadilha: “Meu resíduo é pequeno, mas minha solução está errada”

Para sistemas mal condicionados, você pode ter um resíduo pequeno (|Ax-b|), mas um erro grande na solução (|\hat{x}-x|). Isso não é um bug; é condicionamento.

O que fazer:

• Calcule/estime (\kappa(A)).
• Considere regularização, escalonamento ou reformulação do problema.
• Use métodos baseados em DVS para identificar direções quase nulas.

Armadilha: formar \(A^\top A\) (ou covariância) de forma ingênua

Como observado, isso eleva ao quadrado o condicionamento e pode perder precisão. Além disso, a estimação de covariância pode sofrer cancelamento catastrófico se você a calcular como (E[x^2] - (E[x])^2).

O que fazer:

• Prefira QR/DVS para mínimos quadrados.
• Use algoritmos numericamente estáveis para covariância (centralize os dados primeiro; acumule em maior precisão).

Armadilha: assumir simetria / PSD sem checagens

Uma matriz que “deveria ser” simétrica pode não ser exatamente simétrica devido a erro de ponto flutuante: (A \neq A^\top) no nível de bits.

O que fazer:

• Simetrize explicitamente: (A \leftarrow (A + A^\top)/2) quando apropriado.
• Adicione jitter para matrizes quase PSD.
• Use rotinas do tipo cholesky_ex (quando disponíveis) para detectar falhas de forma robusta.

Armadilha: deficiência de posto e pseudo-inversas

Se (A) é deficiente em posto, muitos problemas se tornam mal postos (múltiplas soluções). Usar uma pseudo-inversa cegamente pode amplificar ruído.

O que fazer:

• Use DVS truncada (descarte valores singulares muito pequenos).
• Use regularização ridge: resolva ((A^\top A + \lambda I)x = A^\top b).

Diagnósticos: como saber o que está dando errado

1) Sempre verifique um resíduo

Para uma solução (Ax=b):

r = b - A @ x_hat
rel_res = np.linalg.norm(r) / (np.linalg.norm(A)*np.linalg.norm(x_hat) + np.linalg.norm(b))

Um resíduo relativo pequeno indica que o algoritmo resolveu bem um problema próximo (erro retroativo pequeno), mas não garante que a solução verdadeira seja precisa se o problema for mal condicionado.

2) Estime o condicionamento (ou pelo menos observe valores singulares)

kappa = np.linalg.cond(A)  # may be expensive for huge matrices

Para cenários em grande escala, aproxime via métodos iterativos ou examine o espectro de (A^\top A) / Fisher / aproximações da Hessiana.

3) Teste de sensibilidade

Perturbe levemente as entradas e veja se as saídas mudam drasticamente:

eps = 1e-7
x1 = np.linalg.lstsq(A, b, rcond=None)[0]
x2 = np.linalg.lstsq(A, b + eps*np.random.randn(*b.shape), rcond=None)[0]
np.linalg.norm(x1 - x2) / np.linalg.norm(x1)

Se uma perturbação minúscula gerar uma mudança enorme, isso é condicionamento, não apenas um “resolvedor com bug”.

Como isso se conecta à prática moderna de AM

• Decomposições estáveis (QR, DVS, Cholesky) importam em AM clássico (modelos lineares, ACP) e dentro de sistemas profundos (branqueamento, blocos de segunda ordem, camadas implícitas).
• Condicionamento explica por que alguns problemas de treinamento são difíceis mesmo quando os gradientes estão corretos — a otimização pode ser lenta ou ruidosa quando a paisagem está mal escalonada (ver Descida do Gradiente).
• Escolhas de precisão (float16/bfloat16 vs float32/float64) interagem com condicionamento: baixa precisão muitas vezes é suficiente para computações bem condicionadas, mas pode falhar de forma espetacular quando (\kappa) é grande.
• Em escala, aproximação (métodos randomizados, resolvedores iterativos) frequentemente é necessária; a álgebra linear numérica fornece ferramentas de erro e estabilidade para tornar essas aproximações confiáveis.

Principais conclusões

• Condicionamento diz respeito ao problema: algumas tarefas são inerentemente sensíveis; (\kappa(A)) grande significa que você deve esperar perda de acurácia.
• Estabilidade diz respeito ao algoritmo: métodos retroativamente estáveis (QR, DVS, LU com pivotamento, Cholesky quando válido) mantêm erros de arredondamento sob controle.
• Evite armadilhas comuns: não compute inversas explícitas, evite equações normais quando a acurácia importa e trate “a matriz é PSD” como uma afirmação numérica que exige validação.
• Para AM em grande escala, o sucesso frequentemente vem da combinação de: formulações estáveis, gestão cuidadosa de precisão e diagnósticos (resíduos + consciência de condicionamento).

Se você for aprofundar este tópico, o próximo passo natural é conectar condicionamento a valores singulares e estrutura de baixa ordem em Autovalores & DVS.

---