Decomposições LU e Cholesky

Fatorar uma matriz é reescrevê-la como produto de fatores mais simples. A ideia central é econômica: fatore uma vez, resolva muitas. Se você precisa resolver $Ax = b$ para vários $b$ diferentes com a mesma $A$, a fatoração transforma um custo $O(n^3)$ por sistema em $O(n^3)$ uma vez e $O(n^2)$ por sistema.

Fatoração LU

$$A = LU$$

com $L$ triangular inferior de diagonal unitária e $U$ triangular superior. É a eliminação gaussiana registrada em forma matricial: $U$ é o resultado da eliminação e $L$ guarda os multiplicadores usados.

Resolver $Ax = b$ vira dois sistemas triangulares:

1. $Ly = b$ por substituição progressiva — $O(n^2)$.
2. $Ux = y$ por substituição regressiva — $O(n^2)$.

Custo da fatoração: $\frac{2}{3}n^3$.

Pivoteamento: PA = LU

A fatoração pura falha quando surge pivô nulo e fica imprecisa com pivô pequeno. Na prática usa-se pivoteamento parcial:

$$PA = LU$$

com $P$ uma matriz de permutação. Isso garante que todos os multiplicadores tenham módulo $\le 1$.

import numpy as np
from scipy.linalg import lu_factor, lu_solve

lu, piv = lu_factor(A)              # fatora uma vez: O(n³)
x1 = lu_solve((lu, piv), b1)        # cada resolução: O(n²)
x2 = lu_solve((lu, piv), b2)

Este padrão — fatorar fora do laço, resolver dentro — é a otimização mais importante de qualquer código que resolve sistemas repetidamente, e é frequentemente esquecida.

Cholesky

Para $A$ simétrica e definida positiva:

$$A = LL^T$$

com $L$ triangular inferior de diagonal positiva.

Vantagens sobre a LU:

• Metade do custo: $\frac{1}{3}n^3$.
• Metade da memória: só um fator.
• Sem pivoteamento: é incondicionalmente estável para matrizes definidas positivas.

import numpy as np
L = np.linalg.cholesky(A)           # falha se A não for definida positiva

O fato de a fatoração falhar quando a matriz não é definida positiva é uma característica útil, não um defeito: é o teste mais barato dessa propriedade. Métodos de otimização usam exatamente isso para detectar que a hessiana deixou de ser definida positiva e que o passo de Newton precisa ser corrigido.

Matrizes simétricas definidas positivas são comuns: matrizes de covariância, matrizes de Gram $A^TA$, hessianas em pontos de mínimo, laplacianos com condição de contorno.

Comparação

Fatoração	Requisito	Custo	Uso
LU	Quadrada	$\frac23 n^3$	Sistemas gerais
Cholesky	Simétrica definida positiva	$\frac13 n^3$	Covariância, mínimos quadrados
QR	$m \ge n$	$\frac43 n^3$	Mínimos quadrados, estabilidade
SVD	Qualquer	$\sim 10n^3$	Posto, pseudoinversa, mal condicionado

A ordem de preferência prática é: use a fatoração mais barata cujas hipóteses a sua matriz satisfaça, exceto quando o condicionamento exigir a mais robusta.

Determinante e log-determinante

Da fatoração sai o determinante de graça:

$$\det(A) = \det(P)^{-1}\prod_i u_{ii}, \qquad \det(A) = \prod_i \ell_{ii}^2 ;\text{(Cholesky)}$$

Para evitar estouro, calcula-se o logaritmo:

L = np.linalg.cholesky(Sigma)
log_det = 2 * np.sum(np.log(np.diag(L)))

Essa expressão aparece diretamente na verossimilhança da distribuição normal multivariada e em processos gaussianos, onde $\Sigma$ pode ter milhares de dimensões.

Amostragem de normal multivariada

Uma aplicação elegante de Cholesky: para amostrar $x \sim \mathcal{N}(\mu, \Sigma)$, basta

$$x = \mu + Lz, \qquad z \sim \mathcal{N}(0, I)$$

pois $\operatorname{Cov}(Lz) = L,I,L^T = LL^T = \Sigma$. Uma fatoração resolve o problema de gerar amostras correlacionadas a partir de ruído independente.

Matrizes esparsas

Para sistemas esparsos grandes, a fatoração sofre preenchimento: zeros viram não zeros durante a eliminação, e o fator pode ser muito mais denso que a matriz original.

Duas respostas:

• Reordenação (mínimo grau aproximado, dissecção aninhada) para reduzir o preenchimento antes de fatorar.
• Métodos iterativos (gradientes conjugados, GMRES), que nunca formam o fator e só precisam do produto matriz-vetor.

Para matrizes muito grandes — milhões de incógnitas em simulação — os métodos iterativos são a única opção viável.

Exemplo trabalhado

$$A = \begin{pmatrix} 4 & 2 \ 2 & 3 \end{pmatrix}$$

Simétrica; testemos Cholesky. Buscamos $L = \begin{pmatrix} \ell_{11} & 0 \ \ell_{21} & \ell_{22}\end{pmatrix}$:

• $\ell_{11}^2 = 4 \Rightarrow \ell_{11} = 2$
• $\ell_{21}\ell_{11} = 2 \Rightarrow \ell_{21} = 1$
• $\ell_{21}^2 + \ell_{22}^2 = 3 \Rightarrow \ell_{22} = \sqrt2$

Como todas as raízes existem e são reais, $A$ é definida positiva. E $\det(A) = (2 \cdot \sqrt2)^2 = 8$, que confere com $4\cdot3 - 2\cdot2 = 8$.

Erros comuns

• Refatorar dentro do laço. O erro de desempenho mais frequente nesta área.
• Cholesky em matriz não definida positiva. Vai falhar — trate a exceção em vez de ignorá-la.
• Ignorar o pivoteamento ao implementar LU manualmente.
• Formar o produto $A^TA$ só para poder usar Cholesky, elevando o condicionamento ao quadrado; prefira QR.
• Fatorar matriz esparsa sem reordenar, causando explosão de memória.

Leituras recomendadas

• Golub e Van Loan, capítulos 3 e 4 — LU e Cholesky em detalhe, incluindo variantes por blocos.
• Trefethen e Bau, lições 20 a 23 — estabilidade da eliminação gaussiana e por que o pivoteamento parcial funciona na prática.
• Davis, Direct Methods for Sparse Linear Systems — a referência sobre fatoração esparsa e reordenação.
• Documentação do SuiteSparse e do scipy.sparse.linalg — as ferramentas de produção.

Gene H. Golub and Charles F. Van Loan (2013). Matrix Computations. Johns Hopkins University Press. DOI: 10.56021/9781421407944. Lloyd N. Trefethen and David Bau III (1997). Numerical Linear Algebra. SIAM. DOI: 10.1137/1.9780898719574. Nicholas J. Higham (2002). Accuracy and Stability of Numerical Algorithms. SIAM. DOI: 10.1137/1.9780898718027.

Referências

1. Gene H. Golub; Charles F. Van Loan (2013). Matrix Computations. Johns Hopkins University Press.
2. Lloyd N. Trefethen; David Bau III (1997). Numerical Linear Algebra. SIAM.
3. Nicholas J. Higham (2002). Accuracy and Stability of Numerical Algorithms. SIAM.