Introducción

La regresión lineal es uno de los algoritmos más básicos y, a la vez, más poderosos del aprendizaje supervisado. A diferencia de muchos tutoriales que se apoyan directamente en scikit‑learn, aquí desglosaremos el algoritmo usando únicamente los pilares de la álgebra lineal (vectores, matrices y operaciones asociadas). De esta forma comprenderás por qué funciona y podrás adaptarlo a casos especiales sin depender de librerías externas.

1. Conceptos clave de Álgebra Lineal

• Vector: \mathbf{x} \in \mathbb{R}^n, una lista ordenada de n valores.
• Matriz: \mathbf{X} \in \mathbb{R}^{m\times n}, colección de m vectores fila (o n vectores columna).
• Producto punto (dot product): \mathbf{a} \cdot \mathbf{b}=\sum_{i=1}^{n} a_i b_i.
• Transpuesta: \mathbf{X}^T intercambia filas y columnas.
• Inversa de una matriz (\mathbf{A}^{-1}) sólo existe si \mathbf{A} es cuadrada y no singular.
• Descomposición en valores singulares (SVD) y pseudoinversa de Moore‑Penrose (\mathbf{X}^{+}) para casos donde \mathbf{X}^T\mathbf{X} no es invertible.

2. Derivación matemática de la regresión lineal

Dados m ejemplos (\mathbf{x}_i, y_i), queremos encontrar un vector de pesos \mathbf{w} que minimice el error cuadrático medio (MSE):

J(\mathbf{w}) = \frac{1}{2m}\| \mathbf{X}\mathbf{w} - \mathbf{y} \|_2^2

Derivando respecto a \mathbf{w} e igualando a cero obtenemos la ecuación normal:

\mathbf{w}^{*} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}

Donde \mathbf{X} incluye una columna de unos para el término independiente (bias).

3. Implementación paso a paso en Python (sin librerías externas)

Esta sección muestra cómo reproducir la ecuación normal usando solo estructuras nativas de Python.

#--- regresion_lineal_puro.py ---
import random
import math
# Función para multiplicar dos matrices (listas de listas)
def mat_mul(A, B):
    result = [[0.0 for _ in range(len(B[0]))] for _ in range(len(A))]
    for i in range(len(A)):
        for j in range(len(B[0])):
            for k in range(len(B)):
                result[i][j] += A[i][k] * B[k][j]
    return result
# Transpuesta de una matriz
def transpose(M):
    return [[M[i][j] for i in range(len(M))] for j in range(len(M[0]))]
# Inversa de una matriz 2x2 (para simplificar la demo)
def inverse_2x2(M):
    a, b = M[0]
    c, d = M[1]
    det = a * d - b * c
    if det == 0:
        raise ValueError('Matriz singular')
    inv_det = 1.0 / det
    return [[ d * inv_det, -b * inv_det],
            [-c * inv_det,  a * inv_det]]
# Generar datos de ejemplo (y = 3 + 2*x + ruido)
random.seed(42)
X_raw = [random.uniform(0, 10) for _ in range(15)]
y = [3 + 2 * x + random.gauss(0, 1) for x in X_raw]
# Construir la matriz de diseño X con columna de bias
X = [[1, x] for x in X_raw]
# Convertir listas a matrices de dimensión (m x 2) y (m x 1)
X_mat = X
y_mat = [[val] for val in y]
# Paso 1: X^T * X
Xt = transpose(X_mat)
XtX = mat_mul(Xt, X_mat)
# Paso 2: (X^T * X)^-1 (solo 2x2 en este caso)
XtX_inv = inverse_2x2(XtX)
# Paso 3: X^T * y
Xty = mat_mul(Xt, y_mat)
# Paso 4: w = (X^T X)^-1 X^T y
w = mat_mul(XtX_inv, Xty)
print('Peso (bias, coef):', w)

Para dimensiones mayores se recomienda usar la pseudoinversa con SVD (implementada en numpy.linalg.pinv), ya que la inversión directa puede fallar por singularidad o por problemas numéricos.

4. Implementación con NumPy y scikit‑learn

import numpy as np
# Datos de ejemplo
X = np.random.rand(100, 1) * 10
y = 3 + 2 * X.squeeze() + np.random.randn(100)
# Añadir bias
X_b = np.c_[np.ones((100, 1)), X]
# Solución cerrada (ecuación normal)
theta_best = np.linalg.inv(X_b.T @ X_b) @ X_b.T @ y
print(theta_best)  # [bias, coef]

from sklearn.linear_model import Ridge
from sklearn.metrics import mean_squared_error, r2_score
model = Ridge(alpha=1.0)  # alpha controla la fuerza de la regularización
model.fit(X, y)
y_pred = model.predict(X)
print('Intercept:', model.intercept_)
print('Coef:', model.coef_)
print('MSE:', mean_squared_error(y, y_pred))
print('R²:', r2_score(y, y_pred))

5. Comparativa de enfoques (dos columnas)

6. Métricas de evaluación y diagnóstico

• Mean Squared Error (MSE): \frac{1}{m}\sum (\hat{y}_i - y_i)^2.
• Root Mean Squared Error (RMSE): \sqrt{MSE} – unidades originales.
• R² (coeficiente de determinación): medida de varianza explicada, 0‑1.
• Residual Plot: inspección visual de la aleatoriedad de los residuos.

Ejemplo de cálculo rápido con NumPy:

residuals = y - X_b @ theta_best
mse = np.mean(residuals**2)
rmse = np.sqrt(mse)
ss_tot = np.sum((y - np.mean(y))**2)
ss_res = np.sum(residuals**2)
r2 = 1 - ss_res/ss_tot
print('MSE:', mse, 'RMSE:', rmse, 'R²:', r2)

7. Buenas prácticas, troubleshooting y seguridad

Escalabilidad

• Para datasets con millones de filas, evita la solución cerrada y usa gradient descent o stochastic gradient descent (SGD) – complejidad O(k·n·p) con k iteraciones.
• Utiliza numpy.linalg.lstsq o scipy.sparse.linalg.lsmr cuando X es dispersa.

Numerical stability

• Escala/normaliza las variables (media 0, varianza 1) para evitar matrices mal condicionadas.
• Prefiere la pseudoinversa (np.linalg.pinv) a la inversa directa.

Diagnóstico de colinealidad

• Calcula el Factor de inflación de la varianza (VIF). Valores >10 indican alta colinealidad.

Seguridad y privacidad

• Si los datos son sensibles, considera differential privacy añadiendo ruido Laplaciano a y antes del entrenamiento.
• Valida que los datos no contengan outliers críticos que puedan sesgar el modelo.

8. Extensiones y tendencias emergentes

La regresión lineal sigue siendo la base de modelos más complejos como Regresión Ridge/Lasso, Generalized Linear Models (GLM) y Deep Learning (capas lineales). En entornos de IA generativa, la interpretación de los pesos lineales se usa para explicar decisiones de modelos de gran escala (XAI).