Descomposición LU: algoritmo, teoría y ejemplos en Python

La descomposición LU (Lower‑Upper) es una de las técnicas más usadas en álgebra lineal numérica para resolver sistemas de ecuaciones, calcular determinantes e invertir matrices de forma eficiente. En este artículo profundizamos en su fundamento matemático, describimos el algoritmo paso a paso y ofrecemos ejemplos de implementación en Python, comparándola con otras factorizaciones como QR y Cholesky.

Fundamentos teóricos

Sea A una matriz cuadrada de dimensión n × n. La descomposición LU busca matrices L (triangular inferior) y U (triangular superior) tal que:

A = L·U

• L tiene 1s en su diagonal (versión Doolittle) o valores libres (versión Crout).
• U es triangular superior sin restricciones en la diagonal.

Si A es singular o necesita pivoting para evitar divisiones por cero, se introduce una matriz de permutación P y la factorización se escribe como P·A = L·U. El pivoting parcial (intercambio de filas) es el más usado por su bajo coste y buena estabilidad numérica.

Algoritmo paso a paso (Doolittle con pivoting parcial)

1. Inicializar P = I_n, L = 0_{n×n} y U = A.
2. Para k = 0 … n‑1:
   • Buscar el pivote p = argmax_{i≥k}|U_{i,k}| y permutar filas k y p en U y P. Si k > 0, aplicar la misma permutación a L (solo columnas <0…k‑1).
   • Para i = k+1 … n‑1:
     • Calcular L_{i,k} = U_{i,k} / U_{k,k}.
     • Actualizar la fila: U_{i, k+1:n} = U_{i, k+1:n} – L_{i,k}·U_{k, k+1:n}.
3. Al final, L tiene 1s en la diagonal y U es triangular superior.

El coste computacional es O(n³), similar a la eliminación de Gauss, pero con la ventaja de que L y U pueden reutilizarse para resolver múltiples sistemas con la misma A.

Implementación en Python

Existen dos caminos habituales:

• Implementación manual con NumPy – útil para comprender cada paso.
• Uso de SciPy – más rápido y robusto, incluye pivoting automático.

1️⃣ Implementación manual (NumPy)

import numpy as np
def lu_decomposition(A):
    n = A.shape[0]
    L = np.eye(n, dtype=float)
    U = A.astype(float).copy()
    P = np.arange(n)
    for k in range(n-1):
        # Pivoting parcial
        piv = np.argmax(np.abs(U[k:, k])) + k
        if piv != k:
            U[[k, piv]] = U[[piv, k]]
            P[[k, piv]] = P[[piv, k]]
            if k > 0:
                L[[k, piv], :k] = L[[piv, 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
# Ejemplo de uso
A = np.array([[2, 3, 1],
              [4, 7, 7],
              [2, 3, 4]], dtype=float)
P, L, U = lu_decomposition(A)
print('P:', P)
print('L:\n', L)
print('U:\n', U)

Este código devuelve el vector de permutación P (no la matriz P completa) y las matrices L y U. Para resolver Ax = b basta a aplicar Pb y luego usar sustitución hacia adelante y hacia atrás.

2️⃣ Uso de SciPy (alta performance)

import numpy as np
from scipy.linalg import lu, lu_solve, lu_factor
A = np.array([[2, 3, 1],
              [4, 7, 7],
              [2, 3, 4]], dtype=float)
P, L, U = lu(A)
print('Matriz P (permutación):\n', P)
print('L:\n', L)
print('U:\n', U)
# Resolver Ax = b con una sola llamada
b = np.array([1, 2, 3], dtype=float)
x = lu_solve(lu_factor(A), b)
print('Solución x:', x)

SciPy maneja internamente el pivoting y devuelve P como matriz de permutación, lo que simplifica la fase de sustitución.

Comparativa con QR y Cholesky

Casos de uso reales

La descomposición LU se emplea en multitud de áreas:

• Ingeniería estructural: resolución de sistemas de ecuaciones estáticas en análisis de elementos finitos.
• Finanzas cuantitativas: cálculo de precios de derivados mediante métodos de Monte Carlo que requieren la solución de sistemas lineales.
• Simulación de circuitos eléctricos: análisis de redes lineales (MNA – Modified Nodal Analysis).
• Machine Learning: algoritmos de regresión lineal y ridge regression cuando se necesita una solución exacta.

Buenas prácticas, optimización y troubleshooting

1. Selección de tipo de LU

Para matrices muy grandes (>10⁴) se recomienda usar versiones basadas en bloques (BLAS/LAPACK) que aprovechan la arquitectura de caché. SciPy ya llama a getrf de LAPACK, que implementa un algoritmo bloqueado.

2. Control de precisión numérica

El factor de crecimiento del pivote (growth factor) indica la pérdida de precisión. Si abs(U).max() es mucho mayor que abs(A).max(), considere usar pivoting completo o scaled partial pivoting.

3. Detección de singularidad

Durante el algoritmo, si abs(U[k,k]) < eps * norm(A), la matriz es prácticamente singular. En SciPy, la excepción LinAlgError indica este caso.

4. Paralelismo

Para entornos de alto rendimiento, emplee numpy.linalg.lu_factor a través de cupy (GPU) o dask.array para distribución en clústeres.

5. Profiling y benchmark

import time, numpy as np, scipy.linalg as la
n = 2000
A = np.random.rand(n, n)
start = time.time()
lu = la.lu_factor(A)
print('Tiempo LU (SciPy):', time.time() - start)

En máquinas con AVX‑512, el factor n³/3 se reduce significativamente.

Conclusión

La descomposición LU sigue siendo una herramienta esencial en el arsenal del ingeniero y científico de datos. Su implementación manual ayuda a comprender los fundamentos de la eliminación gaussiana, mientras que las rutinas optimizadas de SciPy ofrecen velocidad y robustez para producción. Elegir entre LU, QR o Cholesky depende del tipo de matriz y de los requisitos de estabilidad y rendimiento.