1. ¿Qué es la factorización de Cholesky?

La factorización de Cholesky descompone una matriz SPD (Simétrica y Positiva Definida) A ∈ ℝⁿˣⁿ en el producto:

A = L·Lᵀ donde L es una matriz triangular inferior con entradas reales y positivas en la diagonal.

Esta descomposición es la base de muchos algoritmos de optimización, resolución de sistemas lineales y simulaciones Monte Carlo.

2. Fundamento matemático

Para una matriz SPD, la existencia y unicidad de L está garantizada por el teorema de Cholesky. El algoritmo recursivo se define por:

• L₁₁ = sqrt(A₁₁)
• Lᵢ₁ = Aᵢ₁ / L₁₁ (i = 2 … n)
• Para k = 2 … n:
  • Lₖₖ = sqrt(Aₖₖ - Σ_{j=1}^{k-1} Lₖⱼ²)
  • Lᵢₖ = (Aᵢₖ - Σ_{j=1}^{k-1} Lᵢⱼ Lₖⱼ) / Lₖₖ (i = k+1 … n)

El algoritmo tiene complejidad O(n³/3), aproximadamente la mitad del coste de una descomposición LU completa.

3. Implementación básica en Python

Usaremos numpy para operaciones vectorizadas y scipy.linalg para la versión optimizada.

import numpy as np
def cholesky_basic(A):
    """Implementación didáctica del algoritmo de Cholesky (sin pivoteo)."""
    n = A.shape[0]
    L = np.zeros_like(A)
    for k in range(n):
        temp_sum = np.dot(L[k, :k], L[k, :k])
        L[k, k] = np.sqrt(A[k, k] - temp_sum)
        for i in range(k+1, n):
            temp_sum = np.dot(L[i, :k], L[k, :k])
            L[i, k] = (A[i, k] - temp_sum) / L[k, k]
    return L
# Ejemplo de uso
A = np.array([[25, 15, -5],
              [15, 18,  0],
              [-5,  0, 11]], dtype=float)
L = cholesky_basic(A)
print('L =\n', L)
print('Verificación L·Lᵀ =\n', L @ L.T)

Resultado esperado:

L =
[[5.  0.  0.]
 [3.  3.  0.]
 [-1.  1.  3.]]
Verificación L·Lᵀ =
[[25. 15. -5.]
 [15. 18.  0.]
 [-5.  0. 11.]]

4. Versión de producción con scipy.linalg

Para aplicaciones reales, scipy.linalg.cholesky aprovecha BLAS/LAPACK altamente optimizados y ofrece soporte para upper o lower triangular.

from scipy.linalg import cholesky
L = cholesky(A, lower=True)  # devuelve L tal que A = L·Lᵀ
print(L)

Esta función también permite especificar overwrite_a=True para reducir el uso de memoria en matrices muy grandes.

5. Comparativa con otras descomposiciones

6. Buenas prácticas y optimización

• Validar SPD antes de la descomposición: Utiliza np.all(np.linalg.eigvals(A) > 0) o verifica que la diagonal de L sea real y positiva.
• Uso de tipos de dato apropiados: En entornos de alto rendimiento, emplea float32 si la precisión lo permite para reducir la presión en la caché.
• Memoria compartida en procesos paralelos: Cuando se ejecuta en multiprocessing o joblib, comparte la matriz A en memoria de solo lectura para evitar copias.
• Over‑write y alineación: Con scipy.linalg.cholesky(..., overwrite_a=True) y matrices alineadas a 64‑bytes (usando numpy.ndarray(..., align=True)) se consigue un speed‑up del 10‑15% en matrices > 5000×5000.

7. Solución de problemas comunes (troubleshooting)

8. Casos de uso reales

• Gaussian Processes (GP): La inversa de la matriz de covarianza se calcula mediante Cholesky para evitar la costosa inversión directa.
• Optimización convexa (por ejemplo, interior‑point methods): Resolución de sistemas KKT que son SPD.
• Simulación de Monte Carlo: Generación de vectores aleatorios con correlación definida por una matriz de covarianza.
• Control de sistemas: Solución de ecuaciones de Riccati que generan matrices SPD.

9. Escalabilidad y rendimiento en entornos de producción

Para matrices que superan los 10⁵ elementos, se recomiendan:

• Descomposición bloqueada (Block‑Cholesky) disponible en scipy.linalg.blas y bibliotecas como Intel MKL.
• Distribución con Dask o Ray para dividir la matriz en bloques y ejecutar la factorización en clústeres.
• Uso de GPUs con cupy.linalg.cholesky o torch.cholesky para acelerar cálculos en tensores de gran dimensión.

Ejemplo rápido de Block‑Cholesky con Dask:

import dask.array as da
A = da.random.random((5000, 5000), chunks=(1000, 1000))
# Aseguramos SPD
A = A @ A.T + da.eye(5000) * 1e-3
L = da.linalg.cholesky(A)  # ejecuta de forma distribuida