Introducción

El ajuste de curvas es una de las tareas más habituales en ciencia de datos, ingeniería y finanzas. Cuando la relación entre la variable independiente x y la dependiente y se puede describir mediante una combinación lineal de funciones base, el problema se reduce a resolver un sistema de ecuaciones lineales mediante matrices.

En este artículo explicaremos el algoritmo matemático, su implementación paso‑a‑paso en Python usando NumPy y SciPy, y compararemos esta técnica con alternativas como curve_fit y statsmodels. También incluiremos buenas prácticas, consideraciones de rendimiento y un apartado de troubleshooting.

Fundamentos teóricos

Supongamos que queremos ajustar un modelo de la forma:

y ≈ β₀·φ₀(x) + β₁·φ₁(x) + … + βₖ·φₖ(x)

donde φᵢ(x) son funciones base (por ejemplo, potencias de x, funciones trigonométricas, etc.) y βᵢ son los coeficientes que buscamos.

En forma matricial:

Y = X·β + ε

• Y: vector columna de observaciones (n×1).
• X: matriz de diseño (n×(k+1)) cuyas columnas son φᵢ(x).
• β: vector de coeficientes ((k+1)×1).
• ε: error residual.

El criterio de mínimos cuadrados busca β que minimice ‖Y‑Xβ‖₂². La solución analítica es:

β = (XᵀX)⁻¹ Xᵀ Y

Esta ecuación es la base del algoritmo de ajuste por matrices.

Algoritmo paso a paso

1. Definir las funciones base φᵢ(x) y crear la matriz de diseño X.
2. Calcular la transpuesta Xᵀ y el producto XᵀX.
3. Obtener la inversa (o pseudo‑inversa) de XᵀX. En caso de colinealidad, usar np.linalg.pinv.
4. Multiplicar la pseudo‑inversa por XᵀY para obtener β.
5. Evaluar el modelo con Ŷ = X·β y calcular métricas de error (RMSE, R², etc.).

Implementación en Python

Utilizaremos numpy para la manipulación de matrices y matplotlib para visualizar los resultados.

import numpy as np
import matplotlib.pyplot as plt
# Datos de ejemplo (ruido gaussiano)
np.random.seed(42)
x = np.linspace(0, 10, 30)
y_true = 3.5 * x**2 - 2.1 * x + 5.0
noise = np.random.normal(0, 20, size=x.shape)
y = y_true + noise
# 1️⃣ Definir funciones base (polinomio de grado 2)
def design_matrix(x, degree=2):
    """Construye la matriz X con columnas [1, x, x**2, ..., x**degree]"""
    return np.vstack([x**i for i in range(degree + 1)]).T
X = design_matrix(x, degree=2)
# 2️⃣ Solución mediante pseudo‑inversa (más estable que la inversa directa)
beta = np.linalg.pinv(X) @ y
print("Coeficientes estimados:", beta)
# 3️⃣ Predicción y métricas
y_pred = X @ beta
rmse = np.sqrt(np.mean((y - y_pred)**2))
ss_res = np.sum((y - y_pred) ** 2)
ss_tot = np.sum((y - np.mean(y)) ** 2)
r2 = 1 - ss_res / ss_tot
print(f"RMSE = {rmse:.2f}, R² = {r2:.4f}")
# 4️⃣ Visualización
plt.scatter(x, y, label='Datos con ruido', color='steelblue')
plt.plot(x, y_true, label='Modelo verdadero', linewidth=2, color='darkgreen')
plt.plot(x, y_pred, label='Ajuste por matrices', linewidth=2, color='crimson')
plt.legend()
plt.title('Ajuste de curva polinómica usando matrices')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

El script muestra cómo crear la matriz de diseño, obtener los coeficientes mediante la pseudo‑inversa y evaluar el ajuste. En la salida típica obtendremos algo parecido a:

Coeficientes estimados: [ 5.12 -2.08  3.48]
RMSE = 18.73, R² = 0.9621

Los valores están muy cercanos a los parámetros reales (5, -2.1, 3.5).

Comparativa con otras técnicas de ajuste

En la práctica, la elección depende del nivel de complejidad del modelo y de los requisitos de interpretación.

Buenas prácticas, seguridad y rendimiento

1️⃣ Normalización de variables

Escalar x a rango [0, 1] o usar StandardScaler reduce la condición numérica de XᵀX y evita pérdidas de precisión.

2️⃣ Uso de np.linalg.lstsq en vez de la pseudo‑inversa

Para datasets muy grandes (n > 10⁶) lstsq es más rápido y consume menos memoria que pinv porque emplea factorización QR.

3️⃣ Regularización (Ridge)

Si la matriz está mal condicionada, añadir un término de penalización λI (Ridge) estabiliza la solución:

beta = np.linalg.solve(X.T @ X + λ*np.eye(k+1), X.T @ y)

4️⃣ Validación cruzada

Dividir los datos en train/test o usar KFold asegura que el modelo no está sobre‑ajustado.

5️⃣ Seguridad

Cuando los datos provienen de fuentes externas, verifica que no contengan NaN o inf. Un simple filtrado evita errores de LinAlgError: singular matrix:

mask = np.isfinite(x) & np.isfinite(y)
X, y = X[mask], y[mask]

Sección de troubleshooting

• Matrix is singular: Indica colinealidad entre columnas de X. Soluciones: eliminar columnas redundantes, usar pinv o aplicar regularización.
• RMSE muy alto pero R² aceptable: Puede deberse a outliers. Considera usar Huber loss o técnicas robustas.
• Desbordamiento de precisión en float64 con valores extremadamente grandes. Escala los datos o usa np.float128 (si está disponible).
• Lentitud en datasets masivos: Cambia a scipy.sparse.linalg.lsqr con matrices dispersas o emplea dask.array para paralelismo.

Conclusiones

El ajuste de curvas mediante matrices es una herramienta poderosa, especialmente cuando se necesita total control sobre la base funcional y se trabaja en entornos con limitaciones de dependencias. Con NumPy y SciPy podemos obtener soluciones exactas, diagnósticos de condición numérica y extensibilidad a regularizaciones.

Sin embargo, para problemas no lineales complejos o para pipelines de aprendizaje automático con regularización avanzada, vale la pena explorar scikit‑learn, statsmodels o lmfit.

¡Experimenta con distintas bases, valida tus modelos y mantén tus datos limpios! El dominio del ajuste matricial abre la puerta a un modelado estadístico robusto y reproducible.