Phase 01 - Lesson 17

Sistemas Lineales

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

Resolver Ax = b es el problema más antiguo de las matemáticas que aún hace funcionar tu red neuronal.

Tipo: Construir Lenguaje: Python Requisitos previos: Fase 1, Lecciones 01 (Intuición de Álgebra Lineal), 02 (Vectores y Matrices), 03 (Transformaciones de Matrices) Tiempo: ~120 minutos

Objetivos de Aprendizaje

• Resolver Ax = b usando eliminación gaussiana con pivoteo parcial y sustitución hacia atrás
• Factorizar matrices con descomposiciones LU, QR y Cholesky y explicar cuándo cada una es apropiada
• Derivar las ecuaciones normales para mínimos cuadrados y conectarlas con la regresión lineal y ridge
• Diagnosticar sistemas mal condicionados usando el número de condición y aplicar regularización para estabilizarlos

El Problema

Cada vez que entrenas una regresión lineal, resuelves un sistema lineal. Cada vez que calculas un ajuste por mínimos cuadrados, resuelves un sistema lineal. Cada vez que una capa de red neuronal calcula y = Wx + b, está evaluando un lado de un sistema lineal. Cuando agregas regularización, modificas el sistema. Cuando usas procesos gaussianos, factorizas una matriz. Cuando inviertes una matriz de covarianza para la distancia de Mahalanobis, resuelves un sistema lineal.

La ecuación Ax = b aparece en todas partes. A es una matriz de coeficientes conocidos. b es un vector de salidas conocidas. x es el vector de incógnitas que quieres encontrar. En la regresión lineal, A es tu matriz de datos, b es tu vector de objetivos y x es el vector de pesos. El modelo entero se reduce a: encontrar x tal que Ax sea lo más cercano posible a b.

Esta lección construye cada método importante para resolver esa ecuación desde cero. Entenderás por qué algunos métodos son rápidos y otros son estables, por qué algunos funcionan solo para sistemas cuadrados y otros manejan sistemas sobredeterminados, y por qué el número de condición de tu matriz determina si tu respuesta tiene algún significado.

El Concepto

Qué significa Ax = b geométricamente

Un sistema de ecuaciones lineales tiene una interpretación geométrica. Cada ecuación define un hiperplano. La solución es el punto (o conjunto de puntos) donde todos los hiperplanos se intersectan.

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

Pueden ocurrir tres cosas:

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

En forma matricial, "una solución" significa que A es invertible. "Ninguna solución" significa que el sistema es inconsistente. "Infinitas soluciones" significa que A tiene un espacio nulo. La mayoría de los problemas de ML caen en la categoría "ninguna solución exacta" porque tienes más ecuaciones (puntos de datos) que incógnitas (parámetros). Ahí es donde entran los mínimos cuadrados.

Vista por columnas vs vista por filas

Hay dos maneras de leer Ax = b.

Vista por filas. Cada fila de A define una ecuación. Cada ecuación es un hiperplano. La solución es donde todas se intersectan.

Vista por columnas. Cada columna de A es un vector. La pregunta pasa a ser: ¿qué combinación lineal de las columnas de A produce 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.

La vista por columnas es más fundamental. Si b está en el espacio columna de A, el sistema tiene solución. Si no lo está, encuentras el punto más cercano en el espacio columna. Ese punto más cercano es la solución de mínimos cuadrados.

Eliminación gaussiana

La eliminación gaussiana transforma Ax = b en un sistema triangular superior Ux = c que resuelves por sustitución hacia atrás. Es el método más directo.

El 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.

Ejemplo:

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

La eliminación gaussiana cuesta O(n^3) operaciones. Para un sistema 1000x1000, eso es alrededor de mil millones de operaciones de punto flotante. Rápido, pero puedes hacerlo mejor si necesitas resolver múltiples sistemas con la misma A.

Pivoteo parcial: por qué importa

Sin pivoteo, la eliminación gaussiana puede fallar o producir basura. Si un elemento pivote es cero, divides por cero. Si es pequeño, amplificas errores de redondeo.

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.

En la aritmética de punto flotante con precisión limitada, la versión sin pivoteo puede perder dígitos significativos. El pivoteo parcial siempre selecciona el mayor pivote disponible para minimizar la amplificación de errores.

Descomposición LU

La descomposición LU factoriza A en una matriz triangular inferior L y una matriz triangular superior U: A = LU. La matriz L almacena los multiplicadores de la eliminación gaussiana. La matriz U es el resultado de la eliminación.

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 qué factorizar en lugar de solo eliminar? Porque, una vez que tienes L y U, resolver Ax = b para cualquier nuevo b cuesta solo O(n^2):

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

El costo O(n^3) se paga una sola vez durante la factorización. Cada resolución posterior es O(n^2). Si necesitas resolver 1000 sistemas con la misma A pero distintos vectores b, la LU ahorra un factor de 1000/3 en el trabajo total.

Con pivoteo parcial, obtienes PA = LU, donde P es una matriz de permutación que registra los intercambios de filas.

Descomposición QR

La descomposición QR factoriza A en una matriz ortogonal Q y una matriz triangular superior R: A = QR.

Una matriz ortogonal tiene la propiedad Q^T Q = I. Sus columnas son vectores ortonormales. Multiplicar por Q preserva longitudes y á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 es numéricamente más estable que LU para resolver problemas de mínimos cuadrados. El proceso de Gram-Schmidt construye Q columna por columna:

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 paso elimina la componente a lo largo de todos los vectores q anteriores, dejando solo la nueva dirección ortogonal.

Descomposición de Cholesky

Cuando A es simétrica (A = A^T) y definida positiva (todos los valores propios positivos), puedes factorizarla como A = L L^T, donde L es triangular inferior. Esta es la descomposición 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 es dos veces más rápida que LU y requiere la mitad del almacenamiento. Solo funciona para matrices simétricas definidas positivas, pero estas aparecen constantemente:

• Las matrices de covarianza son simétricas semidefinidas positivas (definidas positivas con regularización).
• La matriz de kernel en los procesos gaussianos es simétrica definida positiva.
• La Hessiana de una función convexa en un mínimo es simétrica definida positiva.
• A^T A es siempre simétrica semidefinida positiva.

En los procesos gaussianos, factorizas la matriz de kernel K con Cholesky y luego resuelves K alpha = y para obtener la media predictiva. El factor de Cholesky también proporciona el log-determinante para la verosimilitud marginal: log det(K) = 2 * sum(log(diag(L))).

Mínimos cuadrados: cuando Ax = b no tiene solución exacta

Si A es m x n con m > n (más ecuaciones que incógnitas), el sistema es sobredeterminado. No hay solución exacta. En su lugar, minimizas el error cuadrático:

minimize ||Ax - b||^2

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

El minimizador satisface las ecuaciones normales:

A^T A x = A^T b

Derivación: expande ||Ax - b||^2 = (Ax - b)^T (Ax - b) = x^T A^T A x - 2 x^T A^T b + b^T b. Calcula el gradiente respecto a x e iguálalo a cero: 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.

Ecuaciones normales = regresión lineal

La conexión es exacta. En la regresión lineal, tu matriz de datos X tiene una fila por muestra y una columna por característica. Tu vector de objetivos y tiene una entrada por muestra. El vector de pesos w satisface:

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

Esta es la solución en forma cerrada de la regresión lineal. Cada llamada a sklearn.linear_model.LinearRegression.fit() calcula esto (o un equivalente vía QR o SVD).

Agrega un término de regularización lambda * I a la matriz y obtienes la regresión ridge:

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

La regularización hace que la matriz esté mejor condicionada (más fácil de invertir con precisión) y previene el sobreajuste al encoger los pesos hacia cero. La matriz X^T X + lambda * I siempre es simétrica definida positiva cuando lambda > 0, así que puedes usar Cholesky para resolverla.

Pseudoinversa (Moore-Penrose)

La pseudoinversa A+ generaliza la inversión de matrices a matrices no cuadradas y singulares. Para cualquier matriz A:

x = A+ b

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

Sigma+ se forma tomando el recíproco de cada valor singular no nulo y transponiendo el resultado. Si A = U Sigma V^T, entonces 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

La pseudoinversa proporciona la solución de mínimos cuadrados de norma mínima. Si el sistema tiene:

• Una solución: A+ b la proporciona.
• Ninguna solución: A+ b proporciona la solución de mínimos cuadrados.
• Infinitas soluciones: A+ b proporciona aquella con el menor ||x||.

np.linalg.lstsq y np.linalg.pinv de NumPy usan ambos la SVD internamente.

Número de condición

El número de condición mide cuán sensible es la solución a pequeños cambios en la entrada. Para una matriz A, el número de condición es:

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

donde sigma_max y sigma_min son los valores singulares mayor y menor.

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

Reglas prácticas:

• kappa < 100: seguro, la solución es precisa.
• kappa ~ 10^k: pierdes alrededor de k dígitos de precisión de tu aritmética de punto flotante.
• kappa ~ 10^16 (para float64): la solución no tiene significado. La matriz es efectivamente singular.

En ML, el mal condicionamiento ocurre cuando las características son casi colineales. La regularización (agregar lambda * I) mejora el número de condición de sigma_max / sigma_min a (sigma_max + lambda) / (sigma_min + lambda).

Métodos iterativos: gradiente conjugado

Para sistemas dispersos muy grandes (millones de incógnitas), los métodos directos como LU o Cholesky son demasiado costosos. Los métodos iterativos aproximan la solución mejorando una estimación a lo largo de muchas iteraciones.

El gradiente conjugado (CG) resuelve Ax = b cuando A es simétrica definida positiva. Encuentra la solución exacta en a lo sumo n iteraciones (en aritmética exacta), pero típicamente converge mucho más rápido si los valores propios de A están 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

El CG se usa en:

• Optimización a gran escala (método Newton-CG)
• Resolución de discretizaciones de EDP
• Métodos de kernel donde la matriz de kernel es demasiado grande para factorizar
• Precondicionamiento para otros solucionadores iterativos

La tasa de convergencia depende del número de condición. Los sistemas mejor condicionados convergen más rápido, lo cual es otra razón por la que la regularización ayuda.

El panorama completo: qué método cuándo

Conexión con ML

Cada método de esta lección aparece en ML de producción:

Regresión lineal. La solución en forma cerrada resuelve las ecuaciones normales X^T X w = X^T y. Esto se hace vía Cholesky (si n es pequeño) o QR (si la estabilidad numérica importa) o SVD (si la matriz puede tener rango deficiente).

Regresión ridge. Agrega lambda * I a X^T X. El sistema regularizado (X^T X + lambda * I) w = X^T y siempre se puede resolver vía Cholesky porque X^T X + lambda * I es simétrica definida positiva para lambda > 0.

Procesos gaussianos. La media predictiva requiere resolver K alpha = y, donde K es la matriz de kernel. La factorización de Cholesky de K es el enfoque estándar. La log-verosimilitud marginal usa log det(K) = 2 sum(log(diag(L))).

Inicialización de redes neuronales. La inicialización ortogonal usa la descomposición QR para crear matrices de pesos cuyas columnas son ortonormales. Esto previene el colapso de la señal en redes profundas.

Precondicionamiento. Los optimizadores a gran escala usan Cholesky incompleta o LU incompleta como precondicionadores para los solucionadores de gradiente conjugado.

Ingeniería de características. El número de condición de X^T X indica si tus características son colineales. Si kappa es grande, elimina características o agrega regularización.

Constrúyelo

Paso 1: Eliminación gaussiana con pivoteo 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

Paso 2: Descomposición 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

Paso 3: Descomposición 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

Paso 4: Mínimos cuadrados vía ecuaciones normales

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

Paso 5: Número de condición

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

Úsalo

Juntando las piezas para regresión lineal y regresión ridge con datos reales:

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_}")

Entrégalo

Esta lección produce:

• code/linear_systems.py que contiene implementaciones desde cero de eliminación gaussiana, descomposición LU, descomposición de Cholesky, mínimos cuadrados y regresión ridge
• Una demostración funcional de que las ecuaciones normales y el LinearRegression de sklearn producen los mismos pesos

Ejercicios

1. Resuelve el sistema [[1,2,3],[4,5,6],[7,8,10]] x = [6, 15, 27] usando tu eliminación gaussiana, tu solucionador LU y np.linalg.solve. Verifica que los tres dan la misma respuesta dentro de la tolerancia de punto flotante.

2. Genera una matriz aleatoria 50x5 X y el objetivo y = X @ w_true + ruido. Resuelve para w usando ecuaciones normales, QR (vía np.linalg.qr), SVD (vía np.linalg.svd) y np.linalg.lstsq. Compara las cuatro soluciones. Mide el número de condición de X^T X y explica cómo afecta en qué método confías.

3. Crea una matriz casi singular haciendo dos columnas casi idénticas (por ejemplo, columna 2 = columna 1 + 1e-10 * ruido). Calcula su número de condición. Resuelve Ax = b con y sin regularización (agrega 0.01 * I). Compara las soluciones y los residuos. Explica por qué la regularización ayuda.

4. Implementa el algoritmo del gradiente conjugado para una matriz aleatoria 100x100 simétrica definida positiva. Cuenta cuántas iteraciones le toma converger a la tolerancia 1e-8. Compara con el máximo teórico de n iteraciones.

5. Cronometra tu solucionador de Cholesky vs tu solucionador LU vs np.linalg.solve en matrices simétricas definidas positivas de tamaño 10, 50, 200, 500. Grafica los resultados. Verifica que Cholesky es aproximadamente 2x más rápido que LU.

Términos Clave

Lecturas Adicionales

• MIT 18.06: Linear Algebra (Gilbert Strang) -- el curso definitivo sobre sistemas lineales y factorizaciones de matrices
• Numerical Linear Algebra (Trefethen & Bau) -- la referencia estándar para entender la estabilidad numérica, el condicionamiento y por qué fallan los algoritmos
• Matrix Computations (Golub & Van Loan) -- la referencia enciclopédica para cada algoritmo de matrices
• 3Blue1Brown: Inverse Matrices -- intuición visual para lo que significa resolver Ax = b geométricamente