1. Fundamentos Teóricos

Dados n pares de datos (x_i, y_i), la aproximación por mínimos cuadrados polinomiales busca el polinomio de grado d que minimiza la suma de los residuos al cuadrado:

J(\mathbf{a}) = \sum_{i=1}^{n} \bigl(y_i - p_d(x_i;\mathbf{a})\bigr)^2

donde p_d(x;\mathbf{a}) = a_0 + a_1 x + a_2 x^2 + \dots + a_d x^d y \mathbf{a}= (a_0, a_1, …, a_d) son los coeficientes a determinar.

Al derivar J respecto a cada a_k y igualar a cero se obtiene un sistema lineal de d+1 ecuaciones, conocido como normal equations:

(X^T X)\mathbf{a}=X^T \mathbf{y}

Con X la matriz de Vandermonde:

X = \begin{bmatrix}
    1 & x_1 & x_1^2 & \dots & x_1^d \\
    1 & x_2 & x_2^2 & \dots & x_2^d \\
    \vdots & \vdots & \vdots & & \vdots \\
    1 & x_n & x_n^2 & \dots & x_n^d
\end{bmatrix}

Resolver este sistema nos da los coeficientes óptimos.

2. Paso a Paso del Algoritmo

1. Seleccionar el grado d: balance entre sesgo y varianza.
2. Construir la matriz de Vandermonde X con los valores x_i.
3. Calcular X^T X y X^T y.
4. Resolver el sistema lineal usando factorización QR, SVD o np.linalg.solve para mayor estabilidad numérica.
5. Obtener los coeficientes a_k y construir el polinomio.
6. Evaluar el modelo con métricas como R², RMSE o MAE.

3. Implementación en Python

Usaremos numpy y scipy.linalg. También se muestra la alternativa con numpy.polynomial.Polynomial.fit que internamente emplea np.linalg.lstsq con escalado.

import numpy as np
from scipy.linalg import lstsq
import matplotlib.pyplot as plt
# 1️⃣ Datos de ejemplo (pueden ser reales o sintéticos)
x = np.linspace(-5, 5, 30)
# Función subyacente + ruido
np.random.seed(42)
y = 2.5 - 1.2*x + 0.3*x**2 + np.random.normal(0, 2, size=x.shape)
# 2️⃣ Grado del polinomio (p. ej., 2)
grado = 2
# 3️⃣ Construir la matriz de Vandermonde
X = np.vander(x, N=grado+1, increasing=True)  # columnas: [1, x, x^2, ...]
# 4️⃣ Resolver usando least‑squares (más robusto que la inversa directa)
coef, residuals, rank, s = lstsq(X, y)
print('Coeficientes:', coef)
# 5️⃣ Función que evalúa el polinomio
p = np.poly1d(coef[::-1])  # invertir porque np.poly1d espera orden descendente
# 6️⃣ Visualización
x_fine = np.linspace(x.min(), x.max(), 400)
plt.scatter(x, y, label='Datos', color='steelblue')
plt.plot(x_fine, p(x_fine), 'r--', label=f'Polinomio grado {grado}')
plt.title('Ajuste por Mínimos Cuadrados Polinomial')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.show()

El bloque anterior ilustra todo el flujo, desde la generación de datos hasta la visualización del ajuste.

4. Comparativa con Técnicas Alternativas

5. Buenas Prácticas y Optimización

• Escalado de variables: normalizar x (p. ej., (x - mean)/std) reduce la condición de la matriz de Vandermonde.
• Uso de SVD o QR: np.linalg.lstsq o scipy.linalg.lstsq emplean SVD y son más robustos que la inversa directa.
• Selección del grado: aplicar validación cruzada (k‑fold) para escoger d que maximice R² sin sobre‑ajuste.
• Regularización: añadir un término λ‖a‖² (Ridge) si el grado es alto.
• Vectorización: evitar bucles Python; usar operaciones de numpy para construir X y calcular métricas.

6. Solución de Problemas Comunes

Problema: Matrix is singular or near‑singular

Indicador de que la matriz X^T X está mal condicionada. Soluciones:

• Escalar x antes de crear la Vandermonde.
• Reducir el grado del polinomio.
• Usar np.linalg.lstsq (SVD) que tolera singularidad parcial.

Problema: Sobre‑ajuste evidente (R² ≈ 1 en entrenamiento, bajo en validación)

• Aplicar regularización (Ridge/Lasso).
• Seleccionar un grado menor mediante validación cruzada.
• Considerar splines o modelos no paramétricos.

7. Escalabilidad y Despliegue

Para datasets con millones de observaciones:

• Utiliza scipy.sparse para construir una Vandermonde dispersa cuando el grado es bajo.
• Aplica algoritmos iterativos como Conjugate Gradient sobre X^T X en lugar de factorizaciones directas.
• En entornos de producción, empaqueta la lógica en una función Lambda (AWS) o en un micro‑servicio Docker para responder rápidamente a peticiones de predicción.

8. Conclusión

La aproximación por mínimos cuadrados polinomiales sigue siendo una herramienta esencial para análisis exploratorio y modelado rápido cuando la relación entre variables es razonablemente suave. Con las prácticas de escalado, regularización y uso de algoritmos numéricamente estables, su aplicación es fiable incluso en entornos de producción de alta demanda.