Mínimos cuadrados

La idea

Cuando hay más ecuaciones que incógnitas, como al ajustar una recta a diez puntos medidos, el sistema A·c = y normalmente no tiene solución exacta: los datos traen ruido y ninguna recta pasa por todos. El criterio de mínimos cuadrados renuncia a resolverlo y busca el vector de coeficientes c que hace lo más pequeña posible la suma de los cuadrados de los residuos, es decir de las diferencias entre lo observado y lo predicho. Se elevan al cuadrado para que los errores por exceso y por defecto no se cancelen y para que la función a minimizar sea derivable.

La idea geométrica es una proyección ortogonal. El conjunto de todos los vectores de la forma A·c es el subespacio generado por las columnas de A; el punto de ese subespacio más cercano a y es su proyección, y el residuo que sobra debe ser perpendicular a todas las columnas. Escribir esa condición de perpendicularidad da las ecuaciones normales AᵗA·c = Aᵗ·y, un sistema cuadrado con tantas ecuaciones como parámetros que se resuelve por eliminación de Gauss.

El método no se limita a rectas. Basta elegir bien las columnas de A: con 1 y x se ajusta una recta, añadiendo x² una parábola, y con senos y cosenos un armónico. El modelo sigue siendo lineal en los coeficientes, que es lo único que exige la teoría. Conviene vigilar el grado: aumentar columnas siempre reduce la suma de cuadrados, pero un polinomio demasiado flexible acaba siguiendo el ruido en lugar de la tendencia, y además AᵗA se vuelve cada vez peor condicionada.

Fórmulas clave

• Residuosrᵢ = yᵢ − ŷᵢ, S = Σ rᵢ²
• Ecuaciones normalesAᵗA·c = Aᵗ·y
• Ortogonalidad del residuoAᵗ·(y − A·c) = 0
• Recta de ajustem = Σ(xᵢ − x̄)(yᵢ − ȳ) / Σ(xᵢ − x̄)², b = ȳ − m x̄
• Modelo polinómicoŷ(x) = c₀ + c₁x + c₂x² + … + c_gx^g
• Bondad del ajusteR² = 1 − S / Σ(yᵢ − ȳ)²

Dónde se aplica

En el laboratorio, cualquier curva de calibración —relacionar la absorbancia medida con la concentración de una sustancia— se obtiene por mínimos cuadrados, y su pendiente es el factor de conversión del instrumento. En geodesia y GPS, la posición del receptor se estima ajustando por mínimos cuadrados las distancias redundantes a más satélites de los estrictamente necesarios, lo que reduce el efecto del ruido. En economía, las regresiones que estiman elasticidades son el mismo cálculo, y en aprendizaje automático la regresión lineal con ecuaciones normales es el punto de partida de casi todos los modelos predictivos.

Con programación (Python)

def resolver(A, b):
    n = len(A)
    M = [fila[:] + [b[i]] for i, fila in enumerate(A)]
    for c in range(n):
        p = max(range(c, n), key=lambda r: abs(M[r][c]))
        M[c], M[p] = M[p], M[c]
        for r in range(n):
            if r != c:
                f = M[r][c] / M[c][c]
                for k in range(c, n + 1):
                    M[r][k] -= f * M[c][k]
    return [M[i][n] / M[i][i] for i in range(n)]

def ajuste(xs, ys, grado):
    m = grado + 1
    A = [[x ** j for j in range(m)] for x in xs]
    N = [[sum(A[i][r]*A[i][c] for i in range(len(xs))) for c in range(m)] for r in range(m)]
    t = [sum(A[i][r]*ys[i] for i in range(len(xs))) for r in range(m)]
    return resolver(N, t)

def suma_cuadrados(xs, ys, cs):
    total = 0.0
    for x, y in zip(xs, ys):
        pred = sum(cs[j] * x ** j for j in range(len(cs)))
        total += (y - pred) ** 2
    return total

xs = [0, 1, 2, 3, 4, 5]
ys = [1.1, 2.0, 3.2, 3.9, 5.1, 6.2]
c1 = ajuste(xs, ys, 1)
print("recta   : y =", round(c1[0], 4), "+", round(c1[1], 4), "x   S =", round(suma_cuadrados(xs, ys, c1), 6))
c2 = ajuste(xs, ys, 2)
print("parabola:", [round(v, 5) for v in c2], "  S =", round(suma_cuadrados(xs, ys, c2), 6))
print("suma de residuos (debe ser 0):", round(sum(ys[i] - (c1[0] + c1[1]*xs[i]) for i in range(6)), 9))

La versión en JavaScript, editable y ejecutable en el navegador, está en la página interactiva del tema.

Temas relacionados

Ortogonalidad y Gram-Schmidt · Sistemas lineales y eliminación de Gauss · Regresión lineal simple · Descenso del gradiente