Introducción

La factorización LDLᵀ es una descomposición muy utilizada en análisis numérico y en algoritmos de optimización para matrices simétricas y definidas positivas o semidefinidas. A diferencia de la descomposición de Cholesky (A = LLᵀ), la factorización LDL separa la parte diagonal de la triangular, lo que permite trabajar con matrices que no son estrictamente positivas definidas y evita la necesidad de extraer raíces cuadradas.

Fundamentos matemáticos

Sea A ∈ ℝ^{n×n} una matriz simétrica (A = Aᵀ). La factorización LDL busca matrices:

• L: matriz triangular inferior con unidad en la diagonal (L_{ii}=1).
• D: matriz diagonal (puede contener valores negativos o ceros).

Tal que:

A = L D Lᵀ

El algoritmo procede fila a fila, actualizando L y D mediante operaciones de eliminación gaussiana sin pivoteo (para matrices bien condicionadas).

Algoritmo paso a paso

1. Inicializar L = I_n (identidad) y D = 0_{n×n}.
2. Para k = 0 … n-1:
   • Calcular D_{kk} = A_{kk} - Σ_{s=0}^{k-1} L_{k s}² D_{ss}.
   • Para cada i = k+1 … n-1 actualizar L_{i k} = (A_{i k} - Σ_{s=0}^{k-1} L_{i s} L_{k s} D_{ss}) / D_{kk}.
3. Al final, A = L D Lᵀ.

Si D_{kk}=0 el algoritmo falla; en la práctica se suele aplicar un pivoting parcial o regularizar la matriz añadiendo una pequeña constante a la diagonal.

Comparativa rápida: LDL vs Cholesky

Implementación en Python

Python no incluye una rutina ldl en la API principal de NumPy, pero SciPy ofrece scipy.linalg.ldl a partir de la versión 1.10.0. A continuación se muestra una implementación manual y una alternativa usando SciPy.

Instalación de dependencias

pip install numpy scipy

Implementación manual (educativa)

import numpy as np
def ldl_factorization(A):
    """Devuelve L, D tal que A = L @ D @ L.T.
    A debe ser simétrica.
    """
    n = A.shape[0]
    L = np.eye(n, dtype=A.dtype)
    D = np.zeros((n, n), dtype=A.dtype)
    for k in range(n):
        # Cálculo de D[k,k]
        sum_ld = sum(L[k, s]**2 * D[s, s] for s in range(k))
        D[k, k] = A[k, k] - sum_ld
        if np.isclose(D[k, k], 0):
            raise np.linalg.LinAlgError('Descomposición fallida: pivote cero en D[%d,%d]' % (k, k))
        # Cálculo de los L[i,k]
        for i in range(k+1, n):
            sum_l = sum(L[i, s] * L[k, s] * D[s, s] for s in range(k))
            L[i, k] = (A[i, k] - sum_l) / D[k, k]
    return L, D
# Ejemplo de uso
A = np.array([[4., 12., -16.],
              [12., 37., -43.],
              [-16., -43., 98.]])
L, D = ldl_factorization(A)
print('L =\n', L)
print('D =\n', D)
print('Reconstrucción =\n', L @ D @ L.T)

Uso de SciPy (más rápido y robusto)

from scipy.linalg import ldl
A = np.array([[4., 12., -16.],
              [12., 37., -43.],
              [-16., -43., 98.]])
L, D, perm = ldl(A, lower=True)
print('L =\n', L)
print('D =\n', D)
print('Permutación =', perm)
print('Reconstrucción =\n', L @ D @ L.T)

La función ldl devuelve también un vector de permutación que indica un pivoteo parcial opcional, útil para matrices mal condicionadas.

Ejemplo práctico: Sistema de ecuaciones lineales

Supongamos que queremos resolver Ax = b con la misma A del ejemplo anterior.

import numpy as np
from scipy.linalg import ldl, solve_triangular
A = np.array([[4., 12., -16.],
              [12., 37., -43.],
              [-16., -43., 98.]])
b = np.array([1., 2., 3.])
L, D, perm = ldl(A, lower=True)
# Aplicar permutación a b
b_perm = b[perm]
# Resolución en dos pasos: Ly = b_perm  →  Dz = y  →  Lᵀx = z
y = solve_triangular(L, b_perm, lower=True, unit_diagonal=True)
z = y / np.diag(D)   # D es diagonal, división elemento a elemento
x = solve_triangular(L.T, z, lower=False)
print('Solución x =', x)

El método es comparable en velocidad a la descomposición de Cholesky, pero permite manejar matrices con ceros en la diagonal sin romper la factorización.

Análisis de rendimiento y escalabilidad

• Complejidad temporal: O(n³) en el caso general, idéntica a Cholesky.
• Uso de memoria: Se necesita almacenar L (banda inferior) y D (solo n valores).
• Paralelismo: Operaciones de actualización de fila pueden vectorizarse con BLAS Level‑3; SciPy aprovecha OpenBLAS o MKL bajo el capó.
• Escalabilidad en GPU: Bibliotecas como cuSOLVER ofrecen ldl en CUDA; en Python se puede usar cupy con wrappers personalizados.

Buenas prácticas y troubleshooting

• Verificar simetría: np.allclose(A, A.T). Si la matriz no es exactamente simétrica, symmetrizarla antes de la factorización.
• Condicionamiento: Para matrices mal condicionadas, aplicar regularización (A += ε·I) con ε≈1e-8 para evitar pivotes cero.
• Pivoteo parcial: Utilizar la versión de scipy.linalg.ldl que devuelve la permutación; si perm no es la identidad, reordenar b en sistemas lineales.
• Precisión: En entornos de alta precisión (float64) la factorización es estable; para float32 considere usar dtype=np.float64 en los cálculos intermedios.

Compatibilidad y versiones soportadas

La función scipy.linalg.ldl está disponible a partir de SciPy 1.10.0. Para versiones anteriores, la implementación manual mostrada arriba es completamente compatible con NumPy 1.20+ y Python 3.8+.

En entornos de contenedores (Docker, Podman) se recomienda usar una imagen base python:3.11-slim e instalar las dependencias con:

RUN pip install --no-cache-dir numpy==1.26 scipy==1.12

Conclusión

La factorización LDL es una herramienta esencial cuando se trabaja con matrices simétricas que pueden no ser estrictamente positivas definidas. Su capacidad para exponer la diagonal D la hace valiosa en algoritmos de optimización, análisis estructural y en la solución de sistemas lineales donde la preservación de la escala es crítica. Con las librerías modernas de Python, obtener una descomposición LDL es tan sencillo como una llamada a scipy.linalg.ldl, y la implementación manual sirve como excelente recurso didáctico para comprender los pasos internos del algoritmo.