Análisis de Estabilidad Numérica en Álgebra Lineal
Cómo evaluar y garantizar la precisión de los algoritmos lineales con ejemplos reales en Python.
1. ¿Qué es la estabilidad numérica?
En el contexto del álgebra lineal, la estabilidad numérica mide la capacidad de un algoritmo para producir resultados cercanos al valor exacto a pesar de los errores de redondeo y perturbaciones de datos. Un algoritmo estable no amplifica los errores inherentes al cálculo en punto flotante.
- Estabilidad hacia adelante: el algoritmo devuelve una solución cercana a la solución exacta del problema original.
- Estabilidad hacia atrás: la solución devuelta es la solución exacta de un problema ligeramente perturbado.
En la práctica, la estabilidad está estrechamente relacionada con el número de condición de la matriz y con la forma en que se implementa el algoritmo.
2. Métricas clave para evaluar la estabilidad
Número de condición (cond)
Define cuán sensible es la solución de Ax = b a pequeñas variaciones en A o b. Se calcula como cond(A) = ||A||·||A^{-1}||. Valores grandes (>10⁸) indican problemas potenciales de estabilidad.
import numpy as np
A = np.array([[1, 2], [2, 4.0001]])
cond = np.linalg.cond(A)
print(f"Condición de A: {cond:.2e}")Errores de redondeo y pérdida de significancia
Se cuantifican comparando la solución obtenida con la solución exacta (cuando está disponible) o mediante residuos:
x = np.linalg.solve(A, b)
residual = np.linalg.norm(b - A @ x)
print(f"Residuo: {residual:.2e}")Un residuo pequeño relative al tamaño de b suele indicar buen comportamiento numérico.
3. Algoritmos clásicos y su estabilidad
Métodos Directos
- Eliminación de Gauss (LU): Estable si se aplica pivotación parcial. Sin pivotación, la estabilidad puede colapsar para matrices mal condicionadas.
- Descomposición QR: Muy estable, especialmente útil para sistemas sobre‑determinados y problemas de mínimos cuadrados.
- Descomposición de Cholesky: Requiere matrices simétricas positivas definidas; es numéricamente estable y más eficiente que LU.
Métodos Iterativos
- Gradiente Conjugado (CG): Estable para matrices simétricas positivas definidas; la convergencia depende del número de condición.
- GMRES y BICGSTAB: Apropiados para matrices no simétricas; la estabilidad se mejora con precondicionadores adecuados.
- Jacobi / Gauss‑Seidel: Simples pero sensibles a la condición; útiles como suavizadores en métodos multigrilla.
4. Ejemplos completos en Python
4.1 Solución con LU + Pivotación parcial
import numpy as np
from scipy.linalg import lu_factor, lu_solve
# Matriz mal condicionada
A = np.array([[1e-10, 1], [1, 1]])
b = np.array([1, 2])
# Factorización LU con pivotación
lu, piv = lu_factor(A)
# Resolución del sistema
x = lu_solve((lu, piv), b)
print("Solución LU:", x)
# Evaluación de estabilidad
cond = np.linalg.cond(A)
res = np.linalg.norm(b - A @ x)
print(f"Condición: {cond:.2e}, Residuo: {res:.2e}")
4.2 Descomposición QR para mínimos cuadrados
import numpy as np
from numpy.linalg import qr
# Sistema sobredeterminado (m > n)
A = np.random.rand(8, 3)
b = np.random.rand(8)
Q, R = qr(A, mode='reduced')
# Solución de min‑norm: x = R^{-1} Q^T b
x = np.linalg.solve(R, Q.T @ b)
print("Solución QR (min‑cuadrados):", x)
# Verificación del residuo
res = np.linalg.norm(b - A @ x)
print(f"Residuo: {res:.2e}")
4.3 Gradiente Conjugado con precondicionador Incompleto (IC)
import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import cg, spilu, LinearOperator
# Matriz SPD tridiagonal
n = 1000
k = np.array([np.full(n, -1), np.full(n, 2), np.full(n, -1)])
offset = [-1, 0, 1]
A = diags(k, offset, format='csc')
b = np.ones(n)
# Precondicionador Incompleto LU (IC equivalente para SPD)
M = spilu(A, drop_tol=1e-4)
M_linop = LinearOperator(A.shape, matvec=M.solve)
x, info = cg(A, b, M=M_linop, tol=1e-8)
print('Iteraciones CG:', info)
print('Norma del residuo:', np.linalg.norm(b - A @ x))
5. Buenas prácticas para mantener la estabilidad
- Escalado de matrices: Normaliza filas/columnas para reducir la condición.
- Uso de pivotación en factorizaciones LU y QR.
- Seleccionar el algoritmo adecuado según la estructura (simétrica, positiva definida, esparsa).
- Precondicionadores en métodos iterativos para mejorar la condición efectiva.
- Precisión doble (float64) cuando el número de condición supera 10⁸; considerar precisión extendida (float128) si el hardware lo permite.
- Comprobación de residuos después de cada solución y, si es necesario, aplicar refinamiento iterativo (iterative refinement).
6. Resolución de problemas comunes (troubleshooting)
| Síntoma | Causa probable | Solución recomendada |
|---|---|---|
| Residuo muy grande (>1e-3) | Mala condición o falta de pivotación | Aplicar pivotación parcial; escalar la matriz; usar QR en su lugar. |
| Iteraciones del CG no convergen | Precondicionador inadecuado o matriz no SPD | Verificar simetría; usar GMRES o BICGSTAB con precondicionador ILU. |
| Desbordamiento o NaN en resultados | Overflow por valores extremadamente grandes o bajo flujo | Reescalar datos; utilizar precisión extendida; revisar operaciones de producto punto. |
| Tiempo de ejecución excesivo | Algoritmo O(n³) en matrices grandes | Pasar a métodos iterativos o usar descomposición esparsa (scipy.sparse.linalg). |
7. Compatibilidad, rendimiento y escalabilidad
Los algoritmos presentados están disponibles en NumPy (CPU), SciPy (CPU + sparse) y CuPy (GPU). Para cargas de trabajo a gran escala, considere:
- Uso de BLAS/LAPACK optimizado (OpenBLAS, Intel MKL) para acelerar factorizaciones.
- División de matrices en bloques (Block‑LU, Block‑QR) para aprovechar la memoria caché.
- Implementaciones distribuidas con Dask o PySpark para matrices de varios gigabytes.
Análisis de Estabilidad Numérica en Álgebra Lineal: Algoritmos y Ejemplos Prácticos en Python