Phase 01 - Lesson 17

Sistemas Lineares

This lesson includes a graded coding exercise that runs in your browser, unlocked with lifetime access.

Resolver Ax = b é o problema mais antigo da matemática que ainda faz sua rede neural funcionar.

Tipo: Construir Linguagem: Python Pré-requisitos: Fase 1, Lições 01 (Intuição de Álgebra Linear), 02 (Vetores e Matrizes), 03 (Transformações de Matrizes) Tempo: ~120 minutos

Objetivos de Aprendizagem

• Resolver Ax = b usando eliminação de Gauss com pivoteamento parcial e substituição reversa
• Fatorar matrizes com decomposições LU, QR e Cholesky e explicar quando cada uma é apropriada
• Derivar as equações normais para mínimos quadrados e conectá-las à regressão linear e ridge
• Diagnosticar sistemas mal condicionados usando o número de condição e aplicar regularização para estabilizá-los

O Problema

Toda vez que você treina uma regressão linear, você resolve um sistema linear. Toda vez que você computa um ajuste por mínimos quadrados, você resolve um sistema linear. Toda vez que uma camada de rede neural computa y = Wx + b, ela está avaliando um lado de um sistema linear. Quando você adiciona regularização, você modifica o sistema. Quando você usa processos gaussianos, você fatora uma matriz. Quando você inverte uma matriz de covariância para a distância de Mahalanobis, você resolve um sistema linear.

A equação Ax = b aparece em todo lugar. A é uma matriz de coeficientes conhecidos. b é um vetor de saídas conhecidas. x é o vetor de incógnitas que você quer encontrar. Na regressão linear, A é sua matriz de dados, b é seu vetor de alvos e x é o vetor de pesos. O modelo inteiro se reduz a: encontrar x tal que Ax seja o mais próximo possível de b.

Esta lição constrói cada método importante para resolver essa equação do zero. Você vai entender por que alguns métodos são rápidos e outros são estáveis, por que alguns funcionam apenas para sistemas quadrados e outros lidam com sistemas sobredeterminados, e por que o número de condição da sua matriz determina se a sua resposta tem algum significado.

O Conceito

O que Ax = b significa geometricamente

Um sistema de equações lineares tem uma interpretação geométrica. Cada equação define um hiperplano. A solução é o ponto (ou conjunto de pontos) onde todos os hiperplanos se intersectam.

2x + y = 5          Two lines in 2D.
x - y  = 1          They intersect at x=2, y=1.

graph LR
    A["2x + y = 5"] --- S["Solution: (2, 1)"]
    B["x - y = 1"] --- S

Três coisas podem acontecer:

graph TD
    subgraph "One Solution"
        A1["Lines intersect at a single point"]
    end
    subgraph "No Solution"
        A2["Lines are parallel — no intersection"]
    end
    subgraph "Infinite Solutions"
        A3["Lines are identical — every point is a solution"]
    end

Na forma matricial, "uma solução" significa que A é invertível. "Nenhuma solução" significa que o sistema é inconsistente. "Infinitas soluções" significa que A tem um espaço nulo. A maioria dos problemas de ML cai na categoria "nenhuma solução exata" porque você tem mais equações (pontos de dados) do que incógnitas (parâmetros). É aí que entram os mínimos quadrados.

Visão por colunas vs visão por linhas

Há duas maneiras de ler Ax = b.

Visão por linhas. Cada linha de A define uma equação. Cada equação é um hiperplano. A solução é onde todas elas se intersectam.

Visão por colunas. Cada coluna de A é um vetor. A pergunta passa a ser: qual combinação linear das colunas de A produz b?

A = | 2  1 |    b = | 5 |
    | 1 -1 |        | 1 |

Row picture: solve 2x + y = 5 and x - y = 1 simultaneously.

Column picture: find x1, x2 such that:
  x1 * [2, 1] + x2 * [1, -1] = [5, 1]
  2 * [2, 1] + 1 * [1, -1] = [4+1, 2-1] = [5, 1]   check.

A visão por colunas é mais fundamental. Se b está no espaço-coluna de A, o sistema tem solução. Se não estiver, você encontra o ponto mais próximo no espaço-coluna. Esse ponto mais próximo é a solução de mínimos quadrados.

Eliminação de Gauss

A eliminação de Gauss transforma Ax = b em um sistema triangular superior Ux = c que você resolve por substituição reversa. É o método mais direto.

O algoritmo:

1. For each column k (the pivot column):
   a. Find the largest entry in column k at or below row k (partial pivoting).
   b. Swap that row with row k.
   c. For each row i below k:
      - Compute multiplier m = A[i][k] / A[k][k]
      - Subtract m times row k from row i.
2. Back substitute: solve from the last equation upward.

Exemplo:

Original:
| 2  1  1 | 8 |       R2 = R2 - (2)R1     | 2  1   1 |  8 |
| 4  3  3 |20 |  -->  R3 = R3 - (1)R1 --> | 0  1   1 |  4 |
| 2  3  1 |12 |                            | 0  2   0 |  4 |

                       R3 = R3 - (2)R2     | 2  1   1 |  8 |
                                       --> | 0  1   1 |  4 |
                                           | 0  0  -2 | -4 |

Back substitute:
  -2 * x3 = -4    -->  x3 = 2
  x2 + 2  = 4     -->  x2 = 2
  2*x1 + 2 + 2 = 8 --> x1 = 2

A eliminação de Gauss custa O(n^3) operações. Para um sistema 1000x1000, isso é cerca de um bilhão de operações de ponto flutuante. Rápido, mas você pode fazer melhor se precisar resolver múltiplos sistemas com a mesma A.

Pivoteamento parcial: por que importa

Sem pivoteamento, a eliminação de Gauss pode falhar ou produzir lixo. Se um elemento pivô for zero, você divide por zero. Se ele for pequeno, você amplifica erros de arredondamento.

Bad pivot:                       With partial pivoting:
| 0.001  1 | 1.001 |            Swap rows first:
| 1      1 | 2     |            | 1      1 | 2     |
                                 | 0.001  1 | 1.001 |
m = 1/0.001 = 1000              m = 0.001/1 = 0.001
R2 = R2 - 1000*R1               R2 = R2 - 0.001*R1
| 0.001  1     | 1.001   |      | 1      1     | 2     |
| 0     -999   | -999.0  |      | 0      0.999 | 0.999 |

x2 = 1.000 (correct)            x2 = 1.000 (correct)
x1 = (1.001 - 1)/0.001          x1 = (2 - 1)/1 = 1.000 (correct)
   = 0.001/0.001 = 1.000        Stable because the multiplier is small.

Na aritmética de ponto flutuante com precisão limitada, a versão sem pivoteamento pode perder dígitos significativos. O pivoteamento parcial sempre seleciona o maior pivô disponível para minimizar a amplificação de erros.

Decomposição LU

A decomposição LU fatora A em uma matriz triangular inferior L e uma matriz triangular superior U: A = LU. A matriz L armazena os multiplicadores da eliminação de Gauss. A matriz U é o resultado da eliminação.

A = L @ U

| 2  1  1 |   | 1  0  0 |   | 2  1   1 |
| 4  3  3 | = | 2  1  0 | @ | 0  1   1 |
| 2  3  1 |   | 1  2  1 |   | 0  0  -2 |

Por que fatorar em vez de apenas eliminar? Porque, uma vez que você tem L e U, resolver Ax = b para qualquer novo b custa apenas O(n^2):

Ax = b
LUx = b
Let y = Ux:
  Ly = b    (forward substitution, O(n^2))
  Ux = y    (back substitution, O(n^2))

O custo O(n^3) é pago uma única vez durante a fatoração. Toda resolução subsequente é O(n^2). Se você precisa resolver 1000 sistemas com a mesma A mas vetores b diferentes, a LU economiza um fator de 1000/3 no trabalho total.

Com pivoteamento parcial, você obtém PA = LU, onde P é uma matriz de permutação que registra as trocas de linhas.

Decomposição QR

A decomposição QR fatora A em uma matriz ortogonal Q e uma matriz triangular superior R: A = QR.

Uma matriz ortogonal tem a propriedade Q^T Q = I. Suas colunas são vetores ortonormais. Multiplicar por Q preserva comprimentos e ângulos.

A = Q @ R

Q has orthonormal columns: Q^T Q = I
R is upper triangular

To solve Ax = b:
  QRx = b
  Rx = Q^T b    (just multiply by Q^T, no inversion needed)
  Back substitute to get x.

QR é numericamente mais estável que LU para resolver problemas de mínimos quadrados. O processo de Gram-Schmidt constrói Q coluna por coluna:

Given columns a1, a2, ... of A:

q1 = a1 / ||a1||

q2 = a2 - (a2 . q1) * q1        (subtract projection onto q1)
q2 = q2 / ||q2||                (normalize)

q3 = a3 - (a3 . q1) * q1 - (a3 . q2) * q2
q3 = q3 / ||q3||

R[i][j] = qi . aj    for i <= j

Cada passo remove a componente ao longo de todos os vetores q anteriores, deixando apenas a nova direção ortogonal.

Decomposição de Cholesky

Quando A é simétrica (A = A^T) e positiva definida (todos os autovalores positivos), você pode fatorá-la como A = L L^T, onde L é triangular inferior. Essa é a decomposição de Cholesky.

A = L @ L^T

| 4  2 |   | 2  0 |   | 2  1 |
| 2  5 | = | 1  2 | @ | 0  2 |

L[i][i] = sqrt(A[i][i] - sum(L[i][k]^2 for k < i))
L[i][j] = (A[i][j] - sum(L[i][k]*L[j][k] for k < j)) / L[j][j]    for i > j

Cholesky é duas vezes mais rápida que LU e requer metade do armazenamento. Ela só funciona para matrizes simétricas positivas definidas, mas essas aparecem o tempo todo:

• Matrizes de covariância são simétricas positivas semidefinidas (positivas definidas com regularização).
• A matriz de kernel em processos gaussianos é simétrica positiva definida.
• A Hessiana de uma função convexa em um mínimo é simétrica positiva definida.
• A^T A é sempre simétrica positiva semidefinida.

Em processos gaussianos, você fatora a matriz de kernel K com Cholesky e então resolve K alpha = y para obter a média preditiva. O fator de Cholesky também fornece o log-determinante para a verossimilhança marginal: log det(K) = 2 * sum(log(diag(L))).

Mínimos quadrados: quando Ax = b não tem solução exata

Se A é m x n com m > n (mais equações que incógnitas), o sistema é sobredeterminado. Não há solução exata. Em vez disso, você minimiza o erro quadrático:

minimize ||Ax - b||^2

This is the sum of squared residuals:
  sum((A[i,:] @ x - b[i])^2 for i in range(m))

O minimizador satisfaz as equações normais:

A^T A x = A^T b

Derivação: expanda ||Ax - b||^2 = (Ax - b)^T (Ax - b) = x^T A^T A x - 2 x^T A^T b + b^T b. Calcule o gradiente em relação a x e iguale a zero: 2 A^T A x - 2 A^T b = 0.

Original system (overdetermined, 4 equations, 2 unknowns):
| 1  1 |         | 3 |
| 1  2 | x     = | 5 |       No exact x satisfies all 4 equations.
| 1  3 |         | 6 |
| 1  4 |         | 8 |

Normal equations:
A^T A = | 4  10 |    A^T b = | 22 |
        | 10 30 |            | 63 |

Solve: x = [1.5, 1.7]

This is linear regression. x[0] is the intercept, x[1] is the slope.

Equações normais = regressão linear

A conexão é exata. Na regressão linear, sua matriz de dados X tem uma linha por amostra e uma coluna por feature. Seu vetor de alvos y tem uma entrada por amostra. O vetor de pesos w satisfaz:

X^T X w = X^T y
w = (X^T X)^(-1) X^T y

Esta é a solução de forma fechada da regressão linear. Toda chamada a sklearn.linear_model.LinearRegression.fit() computa isso (ou um equivalente via QR ou SVD).

Adicione um termo de regularização lambda * I à matriz e você obtém a regressão ridge:

(X^T X + lambda * I) w = X^T y
w = (X^T X + lambda * I)^(-1) X^T y

A regularização torna a matriz mais bem condicionada (mais fácil de inverter com precisão) e previne overfitting ao encolher os pesos em direção a zero. A matriz X^T X + lambda * I é sempre simétrica positiva definida quando lambda > 0, então você pode usar Cholesky para resolvê-la.

Pseudoinversa (Moore-Penrose)

A pseudoinversa A+ generaliza a inversão de matrizes para matrizes não quadradas e singulares. Para qualquer matriz A:

x = A+ b

where A+ = V Sigma+ U^T    (computed via SVD)

Sigma+ é formada tomando o recíproco de cada valor singular não nulo e transpondo o resultado. Se A = U Sigma V^T, então A+ = V Sigma+ U^T.

A = U Sigma V^T        (SVD)

Sigma = | 5  0 |       Sigma+ = | 1/5  0  0 |
        | 0  2 |                | 0  1/2  0 |
        | 0  0 |

A+ = V Sigma+ U^T

A pseudoinversa fornece a solução de mínimos quadrados de norma mínima. Se o sistema tem:

• Uma solução: A+ b a fornece.
• Nenhuma solução: A+ b fornece a solução de mínimos quadrados.
• Infinitas soluções: A+ b fornece aquela com o menor ||x||.

np.linalg.lstsq e np.linalg.pinv do NumPy usam ambos a SVD internamente.

Número de condição

O número de condição mede quão sensível é a solução a pequenas mudanças na entrada. Para uma matriz A, o número de condição é:

kappa(A) = ||A|| * ||A^(-1)|| = sigma_max / sigma_min

onde sigma_max e sigma_min são os maiores e menores valores singulares.

Well-conditioned (kappa ~ 1):        Ill-conditioned (kappa ~ 10^15):
Small change in b -->                Small change in b -->
small change in x                    huge change in x

| 2  0 |   kappa = 2/1 = 2          | 1   1          |   kappa ~ 10^15
| 0  1 |   safe to solve            | 1   1+10^(-15) |   solution is garbage

Regras práticas:

• kappa < 100: seguro, a solução é precisa.
• kappa ~ 10^k: você perde cerca de k dígitos de precisão da sua aritmética de ponto flutuante.
• kappa ~ 10^16 (para float64): a solução não tem significado. A matriz é efetivamente singular.

Em ML, o mau condicionamento ocorre quando as features são quase colineares. A regularização (adicionar lambda * I) melhora o número de condição de sigma_max / sigma_min para (sigma_max + lambda) / (sigma_min + lambda).

Métodos iterativos: gradiente conjugado

Para sistemas esparsos muito grandes (milhões de incógnitas), métodos diretos como LU ou Cholesky são caros demais. Métodos iterativos aproximam a solução melhorando um palpite ao longo de muitas iterações.

O gradiente conjugado (CG) resolve Ax = b quando A é simétrica positiva definida. Ele encontra a solução exata em no máximo n iterações (em aritmética exata), mas tipicamente converge muito mais rápido se os autovalores de A estão agrupados.

Algorithm sketch:
  x0 = initial guess (often zero)
  r0 = b - A x0           (residual)
  p0 = r0                 (search direction)

  For k = 0, 1, 2, ...:
    alpha = (rk . rk) / (pk . A pk)
    x_{k+1} = xk + alpha * pk
    r_{k+1} = rk - alpha * A pk
    beta = (r_{k+1} . r_{k+1}) / (rk . rk)
    p_{k+1} = r_{k+1} + beta * pk
    if ||r_{k+1}|| < tolerance: stop

O CG é usado em:

• Otimização de larga escala (método Newton-CG)
• Resolução de discretizações de EDPs
• Métodos de kernel onde a matriz de kernel é grande demais para fatorar
• Pré-condicionamento para outros solvers iterativos

A taxa de convergência depende do número de condição. Sistemas mais bem condicionados convergem mais rápido, o que é mais uma razão pela qual a regularização ajuda.

O panorama completo: qual método quando

Conexão com ML

Cada método desta lição aparece em ML de produção:

Regressão linear. A solução de forma fechada resolve as equações normais X^T X w = X^T y. Isso é feito via Cholesky (se n é pequeno) ou QR (se a estabilidade numérica importa) ou SVD (se a matriz pode ter posto deficiente).

Regressão ridge. Adiciona lambda * I a X^T X. O sistema regularizado (X^T X + lambda * I) w = X^T y é sempre solucionável via Cholesky porque X^T X + lambda * I é simétrica positiva definida para lambda > 0.

Processos gaussianos. A média preditiva requer resolver K alpha = y, onde K é a matriz de kernel. A fatoração de Cholesky de K é a abordagem padrão. A log-verossimilhança marginal usa log det(K) = 2 sum(log(diag(L))).

Inicialização de redes neurais. A inicialização ortogonal usa a decomposição QR para criar matrizes de pesos cujas colunas são ortonormais. Isso previne o colapso de sinal em redes profundas.

Pré-condicionamento. Otimizadores de larga escala usam Cholesky incompleta ou LU incompleta como pré-condicionadores para solvers de gradiente conjugado.

Engenharia de features. O número de condição de X^T X indica se suas features são colineares. Se kappa é grande, remova features ou adicione regularização.

Construa

Passo 1: Eliminação de Gauss com pivoteamento parcial

import numpy as np

def gaussian_elimination(A, b):
    n = len(b)
    Ab = np.hstack([A.astype(float), b.reshape(-1, 1).astype(float)])

    for k in range(n):
        max_row = k + np.argmax(np.abs(Ab[k:, k]))
        Ab[[k, max_row]] = Ab[[max_row, k]]

        if abs(Ab[k, k]) < 1e-12:
            raise ValueError(f"Matrix is singular or nearly singular at pivot {k}")

        for i in range(k + 1, n):
            m = Ab[i, k] / Ab[k, k]
            Ab[i, k:] -= m * Ab[k, k:]

    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (Ab[i, -1] - Ab[i, i+1:n] @ x[i+1:n]) / Ab[i, i]

    return x

Passo 2: Decomposição LU

def lu_decompose(A):
    n = A.shape[0]
    L = np.eye(n)
    U = A.astype(float).copy()
    P = np.eye(n)

    for k in range(n):
        max_row = k + np.argmax(np.abs(U[k:, k]))
        if max_row != k:
            U[[k, max_row]] = U[[max_row, k]]
            P[[k, max_row]] = P[[max_row, k]]
            if k > 0:
                L[[k, max_row], :k] = L[[max_row, k], :k]

        for i in range(k + 1, n):
            L[i, k] = U[i, k] / U[k, k]
            U[i, k:] -= L[i, k] * U[k, k:]

    return P, L, U

def lu_solve(P, L, U, b):
    n = len(b)
    Pb = P @ b.astype(float)

    y = np.zeros(n)
    for i in range(n):
        y[i] = Pb[i] - L[i, :i] @ y[:i]

    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (y[i] - U[i, i+1:] @ x[i+1:]) / U[i, i]

    return x

Passo 3: Decomposição de Cholesky

def cholesky(A):
    n = A.shape[0]
    L = np.zeros_like(A, dtype=float)

    for i in range(n):
        for j in range(i + 1):
            s = A[i, j] - L[i, :j] @ L[j, :j]
            if i == j:
                if s <= 0:
                    raise ValueError("Matrix is not positive definite")
                L[i, j] = np.sqrt(s)
            else:
                L[i, j] = s / L[j, j]

    return L

Passo 4: Mínimos quadrados via equações normais

def least_squares_normal(A, b):
    AtA = A.T @ A
    Atb = A.T @ b
    return gaussian_elimination(AtA, Atb)

def ridge_regression(A, b, lam):
    n = A.shape[1]
    AtA = A.T @ A + lam * np.eye(n)
    Atb = A.T @ b
    L = cholesky(AtA)
    y = np.zeros(n)
    for i in range(n):
        y[i] = (Atb[i] - L[i, :i] @ y[:i]) / L[i, i]
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (y[i] - L.T[i, i+1:] @ x[i+1:]) / L.T[i, i]
    return x

Passo 5: Número de condição

def condition_number(A):
    U, S, Vt = np.linalg.svd(A)
    return S[0] / S[-1]

Use

Juntando as peças para regressão linear e regressão ridge em dados reais:

np.random.seed(42)
X_raw = np.random.randn(100, 3)
w_true = np.array([2.0, -1.0, 0.5])
y = X_raw @ w_true + np.random.randn(100) * 0.1

X = np.column_stack([np.ones(100), X_raw])

w_ols = least_squares_normal(X, y)
print(f"OLS weights (ours):    {w_ols}")

w_np = np.linalg.lstsq(X, y, rcond=None)[0]
print(f"OLS weights (numpy):   {w_np}")
print(f"Max difference: {np.max(np.abs(w_ols - w_np)):.2e}")

w_ridge = ridge_regression(X, y, lam=1.0)
print(f"Ridge weights (ours):  {w_ridge}")

from sklearn.linear_model import Ridge
ridge_sk = Ridge(alpha=1.0, fit_intercept=False)
ridge_sk.fit(X, y)
print(f"Ridge weights (sklearn): {ridge_sk.coef_}")

Entregue

Esta lição produz:

• code/linear_systems.py contendo implementações do zero de eliminação de Gauss, decomposição LU, decomposição de Cholesky, mínimos quadrados e regressão ridge
• Uma demonstração funcional de que as equações normais e o LinearRegression do sklearn produzem os mesmos pesos

Exercícios

1. Resolva o sistema [[1,2,3],[4,5,6],[7,8,10]] x = [6, 15, 27] usando sua eliminação de Gauss, seu solver LU e np.linalg.solve. Verifique que os três fornecem a mesma resposta dentro da tolerância de ponto flutuante.

2. Gere uma matriz aleatória 50x5 X e o alvo y = X @ w_true + ruído. Resolva para w usando equações normais, QR (via np.linalg.qr), SVD (via np.linalg.svd) e np.linalg.lstsq. Compare as quatro soluções. Meça o número de condição de X^T X e explique como ele afeta em qual método você confia.

3. Crie uma matriz quase singular tornando duas colunas quase idênticas (por exemplo, coluna 2 = coluna 1 + 1e-10 * ruído). Compute seu número de condição. Resolva Ax = b com e sem regularização (adicione 0.01 * I). Compare as soluções e os resíduos. Explique por que a regularização ajuda.

4. Implemente o algoritmo do gradiente conjugado para uma matriz aleatória 100x100 simétrica positiva definida. Conte quantas iterações ele leva para convergir à tolerância 1e-8. Compare com o máximo teórico de n iterações.

5. Cronometre seu solver de Cholesky vs seu solver LU vs np.linalg.solve em matrizes simétricas positivas definidas de tamanho 10, 50, 200, 500. Plote os resultados. Verifique que Cholesky é cerca de 2x mais rápido que LU.

Termos-Chave

Leitura Complementar

• MIT 18.06: Linear Algebra (Gilbert Strang) -- o curso definitivo sobre sistemas lineares e fatorações de matrizes
• Numerical Linear Algebra (Trefethen & Bau) -- a referência padrão para entender estabilidade numérica, condicionamento e por que algoritmos falham
• Matrix Computations (Golub & Van Loan) -- a referência enciclopédica para cada algoritmo de matrizes
• 3Blue1Brown: Inverse Matrices -- intuição visual para o que significa resolver Ax = b geometricamente