CO.DE Academy — Aprende tecnología con práctica real

CO.DE Academy — Aprende tecnología con práctica real

Cuando el sistema no tiene solución (y eso está bien)

En cualquier problema con datos medidos hay más ecuaciones que incógnitas: 200 observaciones, 3 parámetros. El sistema $A\mathbf{c} = \mathbf{b}$ exige que $\mathbf{b}$ caiga exactamente en $\operatorname{Col}(A)$ , un subespacio de dimensión 3 dentro de $\mathbb{R}^{200}$ . Con ruido de medida, eso no ocurre nunca.

La reacción correcta no es rendirse ni «forzar» el sistema, sino cambiar la pregunta:

$\text{en vez de } A\mathbf{c} = \mathbf{b}, \qquad \min_{\mathbf{c}} \lVert \mathbf{b} - A\mathbf{c}\rVert^2$

Geométricamente: buscar el punto de $\operatorname{Col}(A)$ más cercano a $\mathbf{b}$ . Y ya sabemos, desde la lección 3, cómo se encuentra el punto más cercano de un subespacio: proyectando.

Proyección sobre un subespacio

La lección 3 proyectaba sobre una recta imponiendo que el residuo fuera perpendicular a un vector. Ahora el subespacio tiene dimensión $n$ y hay que imponer que el residuo sea perpendicular a todas las columnas a la vez. Eso es exactamente $A^{\mathsf{T}}\mathbf{e} = \mathbf{0}$ :

$A^{\mathsf{T}}(\mathbf{b} - A\hat{\mathbf{c}}) = \mathbf{0} \quad\Longleftrightarrow\quad A^{\mathsf{T}}A\,\hat{\mathbf{c}} = A^{\mathsf{T}}\mathbf{b}$

Estas son las ecuaciones normales. Obsérvese lo que ha pasado: un sistema de 200 ecuaciones incompatible se ha convertido en uno de 3 ecuaciones con 3 incógnitas y solución única. Proyectar reduce el problema a la dimensión del modelo.

Si las columnas de $A$ son independientes, $A^{\mathsf{T}}A$ es invertible —se demuestra viendo que $A^{\mathsf{T}}A\mathbf{x} = \mathbf{0}$ implica $\lVert A\mathbf{x}\rVert^2 = 0$ , luego $A\mathbf{x} = \mathbf{0}$ , luego $\mathbf{x} = \mathbf{0}$ — y entonces

$\hat{\mathbf{c}} = (A^{\mathsf{T}}A)^{-1}A^{\mathsf{T}}\mathbf{b}, \qquad \mathbf{p} = \underbrace{A(A^{\mathsf{T}}A)^{-1}A^{\mathsf{T}}}_{P}\,\mathbf{b}$

La matriz $P$ es la generalización directa de la de la lección 3, y cumple lo mismo: $P^{\mathsf{T}} = P$ y $P^2 = P$ .

1 de 6

Por qué el mínimo es único y por qué no hay que buscarlo

La función que se minimiza es

$f(\mathbf{c}) = \lVert\mathbf{b} - A\mathbf{c}\rVert^2 = \mathbf{c}^{\mathsf{T}}(A^{\mathsf{T}}A)\mathbf{c} - 2\mathbf{b}^{\mathsf{T}}A\mathbf{c} + \mathbf{b}^{\mathsf{T}}\mathbf{b}$

una forma cuadrática (lección 11) con matriz $A^{\mathsf{T}}A$ definida positiva cuando las columnas son independientes. Su gráfica es un cuenco convexo: no hay mínimos locales, ni mesetas, ni sillas. Por eso las curvas de nivel del laboratorio son elipses concéntricas alrededor de un único punto, y por eso el ajuste lineal es tan estable comparado con casi cualquier ajuste no lineal: aquí no hay que buscar, hay que resolver.

La forma de esas elipses no es un detalle estético. Cuanto más alargadas, peor condicionado está el problema y más sensible es la solución al ruido — el número de condición de la lección 7, ahora visible.

Bases ortonormales: por qué se persiguen

Si las columnas de $A$ fueran ortonormales —llamémoslas $Q$ —, entonces $Q^{\mathsf{T}}Q = I$ y todo se simplifica de golpe:

$\hat{\mathbf{c}} = Q^{\mathsf{T}}\mathbf{b}, \qquad P = QQ^{\mathsf{T}}$

Sin inversas, sin sistemas: cada coordenada es un producto punto, y cambiar una no afecta a las demás. Ese desacoplamiento es lo que hace útiles a las bases ortonormales en todas partes: los senos y cosenos de Fourier en Señales y Sistemas, las componentes principales en la lección 12, las bases de wavelets en compresión de imagen.

Gram-Schmidt: fabricar una base ortonormal

Dada una base cualquiera $\{\mathbf{a}_1,\ldots,\mathbf{a}_n\}$ , el procedimiento construye una ortonormal restando a cada vector sus proyecciones sobre los anteriores —es decir, quitándole todo lo que ya estaba representado:

$\mathbf{v}_k = \mathbf{a}_k - \sum_{i<k} (\mathbf{q}_i^{\mathsf{T}}\mathbf{a}_k)\,\mathbf{q}_i, \qquad \mathbf{q}_k = \frac{\mathbf{v}_k}{\lVert\mathbf{v}_k\rVert}$

Si en algún paso $\mathbf{v}_k = \mathbf{0}$ , ese $\mathbf{a}_k$ era combinación de los anteriores: el procedimiento detecta la dependencia lineal sin buscarla.

Guardando los coeficientes usados se obtiene, gratis, la factorización QR:

$A = QR, \qquad R = Q^{\mathsf{T}}A \text{ triangular superior}$

Y con ella los mínimos cuadrados se resuelven sin formar nunca $A^{\mathsf{T}}A$ :

$A^{\mathsf{T}}A\hat{\mathbf{c}} = A^{\mathsf{T}}\mathbf{b} \Longrightarrow R^{\mathsf{T}}Q^{\mathsf{T}}QR\,\hat{\mathbf{c}} = R^{\mathsf{T}}Q^{\mathsf{T}}\mathbf{b} \Longrightarrow R\,\hat{\mathbf{c}} = Q^{\mathsf{T}}\mathbf{b}$

un sistema triangular, resoluble por sustitución hacia atrás. Esto importa mucho más de lo que parece: formar $A^{\mathsf{T}}A$ eleva al cuadrado el número de condición, $\kappa(A^{\mathsf{T}}A) = \kappa(A)^2$ . Un problema con $\kappa(A) = 10^8$ —perfectamente resoluble por QR— se vuelve irresoluble en doble precisión si se pasa por las ecuaciones normales. Por eso las bibliotecas serias resuelven mínimos cuadrados con QR o con SVD, nunca con las ecuaciones normales.

Ejemplo resuelto: ajustar una recta a mano

Ajustemos $y = c_0 + c_1x$ a los puntos $(1,2)$ , $(2,3)$ y $(3,5)$ . La matriz de diseño y el término independiente son

$A = \begin{bmatrix}1&1\\1&2\\1&3\end{bmatrix}, \qquad \mathbf{b} = \begin{bmatrix}2\\3\\5\end{bmatrix}$

Las ecuaciones normales:

$A^{\mathsf{T}}A = \begin{bmatrix}3&6\\6&14\end{bmatrix}, \qquad A^{\mathsf{T}}\mathbf{b} = \begin{bmatrix}10\\23\end{bmatrix}$

Resolviendo: $3c_0 + 6c_1 = 10$ y $6c_0 + 14c_1 = 23$ . Restando el doble de la primera a la segunda: $2c_1 = 3$ , luego $c_1 = 1.5$ y $c_0 = (10 - 9)/3 = 0.3333$ . La recta ajustada es $y = 0.3333 + 1.5x$ .

Los valores predichos son $1.8333$ , $3.3333$ y $4.8333$ , y los residuos $0.1667$ , $-0.3333$ , $0.1667$ . Dos comprobaciones que siempre deben cumplirse: los residuos suman cero (porque la primera columna de $A$ es de unos, y el residuo es ortogonal a ella) y su producto punto con la columna $x$ también es cero: $0.1667(1) - 0.3333(2) + 0.1667(3) = 0$ .

Un problema de mínimos cuadrados tiene $\kappa(A) \approx 10^{8}$ . ¿Qué método conviene usar?

Lo que mínimos cuadrados promete y lo que no

Minimizar la suma de cuadrados tiene una justificación estadística sólida —es el estimador de máxima verosimilitud si el ruido es gaussiano e independiente— y una debilidad conocida: castiga los errores grandes de forma cuadrática, así que un solo dato atípico puede arrastrar toda la recta. Compruébalo en el laboratorio activando el outlier.

Dos avisos más, que valen para cualquier ajuste:

Añadir parámetros siempre baja el RSS. Pasar de recta a cúbica mejora el ajuste por construcción, aunque el modelo esté empezando a ajustar el ruido. Comparar modelos exige penalizar la complejidad, no mirar solo el error.
Lineal se refiere a los parámetros, no a la curva. Ajustar una parábola $c_0 + c_1x + c_2x^2$ es un problema lineal: las incógnitas entran linealmente aunque $x$ aparezca al cuadrado. Ajustar $c_0e^{c_1x}$ no lo es.

Playground · python

Lo que este modelo deja fuera

El Gram-Schmidt clásico que hemos escrito es el correcto en aritmética exacta y inestable en coma flotante: los vectores dejan de ser ortogonales según avanza. Las implementaciones reales usan la variante modificada o, mejor, reflexiones de Householder. Y si las columnas de $A$ son dependientes, $A^{\mathsf{T}}A$ no es invertible y no hay solución única; la elección canónica —la de norma mínima— exige la pseudoinversa de la lección 12. En el laboratorio, los datos vienen de un generador de semilla fija con un modelo cuadrático conocido, de modo que puedes comparar los coeficientes ajustados con los verdaderos.

Hasta aquí hemos estudiado qué hace una matriz a los subespacios. En la lección 10 preguntamos otra cosa: ¿hay direcciones que la matriz respete, que solo estire sin girar?