Métodos potenciales de prospección, FCAG, 2023.
Aquí haremos una revisión de distintas alternativas de resolver el problema de ajustar una recta a una serie de puntos. Las resultados que se desprenden de esta actividad tienen una importancia general en los métodos de inversión. Especialmente en la temática del aprendizaje automático (*machine learning*), que actualmente ocupa un lugar destacado en la geofísica y en distintas ramas de las ciencia y tecnología.
$\newcommand{\mean}[1]{<\!{#1}\!>}$ $\newcommand{\B}[1]{\mathbf{#1}}$ Generamos el dato que usaremos a lo largo de toda la actividad:
a0, b0 = 10., 5.
y_dato = a0 + b0 * x + ruido
Hemos generado las "observaciones" ${\mathbf{y}}= (y_1,\cdots,y_k,\cdots,y_n)^T$ en los puntos dados en ${\mathbf{x}}= (x_1,\cdots,x_k,\cdots,x_n)^T$.
El nombre regresión lineal sugiere que lo que haremos es devolver las observaciones a la variedad (manifold) donde pertenecen, en este caso una recta. Las observaciones se han alejado de su variedad debido al ruido aleatorio agregado a ellas.
$\def\mean#1{\left< #1 \right>}$ A modo de introducción, consideramos el caso donde queremos ajustar a las observaciones $\mathbf{y}$ un modelo dado por un valor constante $a$. Entonces cada observación $y_k$ puede aproximarse por \begin{equation} \hat{y}_k = a, \end{equation} el modelo es independiente de la posición $x_k$ donde se obtuvo $y_k$.
Para estimar un valor de $a$, minimizaremos la siguiente función objetivo:
\begin{equation} q(a) := \dfrac{1}{2}|| \mathbf{y} - \hat{\mathbf{y}} ||^2, \end{equation}con $||\mathbf{y}||=\sqrt{\sum_{j} y_j^2}$ la norma Euclidea o $L_2$.
Calculando $\dfrac{\mathrm{d}}{\mathrm{d}a}q$ e igualando el resultado a $0$, resulta:
\begin{equation} \hat{a} = \frac{1}{N}\sum_{k=1}^{N}y_{k} = \mean{\mathbf{y}}, \end{equation}donde $\mean{\mathbf{y}}$ es el promedio de los $N$ puntos $y_k$.
Minimizar la función de costo en esta situación nos condujo al operador de promedio simple como la solución que minimiza el error entre los datos y el modelo propuesto.
a = 12.414 (redondeo a tres decimales)
Calculemos la energía del residuo entre las observaciones $\B{y}$, y el modelo propuesto $\mean{\B{y}}$ :
def SSR(u,v):
# Sum of Square Residuals (SSR)
return np.sum(np.square(u-v))
donde $\textrm{SSR}(\B{u},\B{v})$ es la suma de los cuadrados de la diferencias entre los elementos de $\B{u}$ y $\B{v}$.
SSR(y,<y>) = 298.56
En lo que sigue usaremos el vector columna $\mathbf{m}$ para listar los parámetros. En este caso es $\mathbf{m}=(a)$. También denotaremos ${\hat{\mathbf{y}}(\mathbf{m)}}$ como nuestro modelo que depende de los parámetros en $\mathbf{m}$.
$\def\mean#1{\left< #1 \right>}$ El modelo propuesto es ahora:
\begin{equation} \hat{y}_k = a + bx_k. \end{equation}Notemos que buscamos ajustar $M=2$ parámetros (pendiente y ordenada al orígen) a $N$ observaciones. Dado que $N>M$ tenemos un sistema de ecuaciones sobredeterminado.
La función objetivo es entonces función de dos parámetros,
\begin{equation} q = q(\mathbf{m}) = q(a,b) := \dfrac{1}{2}|| \mathbf{y} - \hat{\mathbf{y}} ||^2 = \frac{1}{2}\sum_k (y_k-\hat{y}_k)^2 . \end{equation}Los valores $\hat{a},\hat{b}$ que minimizan $q$ los obtenemos calculando las derivadas parciales de $q$ respecto de $a$ y $b$ e igualando a cero. En este caso, existen expresiones analíticas para los parámetros del ajuste de una recta:
\begin{align} \hat{a} &=\mean{\mathbf{y}}-\mean{\mathbf{x}}\hat{b},\\ \hat{b}&=\frac{\mean{\mathbf{x}\circ\mathbf{y}}-\mean{\mathbf{x}}\mean{\mathbf{y}}}{\mean{\mathbf{x}\circ\mathbf{x}}-\mean{\mathbf{x}}^2}. \nonumber \end{align}Los vectores $\mathbf{x}\circ\mathbf{y}$ y $\mathbf{x}\circ\mathbf{x}$ son productos de Hadamard (productos punto a punto), que resultan en vectores de elementos $x_{k}y_{k}$ y $x_{k}^2$, respectivamente. La notación $\mean{\cdot}$ indica tomar el promedio.
En este caso el vector de parámetros es $\mathbf{m}=(a,b)^T$. Notemos que hemos calculado $\dfrac{\partial}{\partial a}q$ y $\dfrac{\partial}{\partial b}q$ e igualado a cero ambas expresiones para obtener $a$ y $b$. Es decir, hemos resuelto el sistema de ecuaciones lineales sobredeterminado:
\begin{equation} (\dfrac{\partial}{\partial a},\dfrac{\partial}{\partial b})\,q = {\nabla}_{\mathbf{m}}{q}= \mathbf{0}. \end{equation}Para ello, hemos calculamos el gradiente de la función de costo respecto de los parámetros del modelo ${\nabla}_{\mathbf{m}}{q}$.
Programamos la solución analítica de los coeficientes de la regresión lineal, calculamos y graficamos:
a = 9.874 (redondeo a tres decimales) b = 5.277 (redondeo a tres decimales)
Calculemos la energía del residuo entre las observaciones $\B{y}$, y el modelo propuesto $\B{a}+b\B{x}$ :
SSR(y,a+b*x) = 95.4
Podemos calcular el estadístico $R^2$ entre el modelo lineal obtenido, $\mathbf{a}+b\mathbf{x}$, y el modelo que consiste en considerar el valor medio de las mediciones $\mean{\B{y}}$:
$$R^2 = \dfrac{\textrm{SSR}(\mathbf{y},\mean{\B{y}})-\textrm{SSR}(\B{y},\B{a}+b\B{x})}{\textrm{SSR}(\B{y},\mean{\B{y}})},$$El ajuste resulta en un R2 = 0.68
Interpretamos este resultado de la siguiente forma: suponer un modelo lineal representa una reducción del 68% del error de ajuste de suponer el valor medio como modelo de las observaciones.
Denotamos el modelo lineal como
\begin{equation} \hat{\mathbf{y}}(\mathbf{x};\mathbf{m}):=A(\mathbf{x})\,\mathbf{m}=A\mathbf{m}, \end{equation}
donde $A$ ($N\times 2$) se conoce como matriz de diseño. La matriz $A$ tiene el valor $1$ en la primer columna y a los valores de $x_k$ en la segunda columna:
\begin{equation} A := \left[ \mathbf{1}\quad \mathbf{x} \right], \end{equation}
es decir $A_{i,j}=x_i^{j-1}$ para $i=1,\cdots,N$ y $j=1,2$.
Construimos la función objetivo
\begin{equation} q(\mathbf{m}) := \dfrac{1}{2}|| \mathbf{y} - \hat{\mathbf{y}} ||^2=\dfrac{1}{2}|| \mathbf{y} - A\mathbf{m} ||^2, \end{equation}
y resolvemos algebraicamente $\nabla_\mathbf{m}q=\mathbf{0}$ para $\mathbf{m}$. Del paso anterior resulta la conocida solución al sistema normal de ecuaciones:
\begin{equation} \mathbf{m} = (A^T A)^{-1} A^T\mathbf{y}, \end{equation}
siempre y cuando $(A^T A)^{-1}$ exista. En este caso como la matriz $A^TA$ es $2\times 2$, tenemos una forma analítica de calcular $A^{-1}$. ¿La recuerdan o la pueden deducir?
Programamos las subrutinas necesarias,
y resolvemos las ecuaciones normales,
a = 9.874 (redondeo a tres decimales) b = 5.277 (redondeo a tres decimales)
A modo de comparación, utilizamos subrutinas especializadas de SciPy para resolver el sistema de ecuaciones normales:
from scipy.stats import linregress
sol = linregress(x,y)
print("a =", sol.intercept)
print("b =", sol.slope)
a = 9.874 (redondeo a tres decimales) b = 5.277 (redondeo a tres decimales)
La función objetivo
\begin{equation} q := \dfrac{1}{2}|| \mathbf{y} - \hat{\mathbf{y}} ||^2 = \dfrac{1}{2}|| \mathbf{y} - A{\mathbf{m}} ||^2, \end{equation}tiene por gradiente:
\begin{equation} {\nabla}_{\mathbf{m}}{q}= {A^T}({A}\mathbf{m}-\mathbf{y}). \end{equation}El método de descenso del gradiente (gradient descent) con tamaño de paso fijo $\epsilon$ que implementamos es:
Si utilizamos un valor inicial y $\epsilon$ adecuados, utilizando este método recurrimos a un elevado número de iteraciones para obtener una solución aceptable. En oposición a resolver el sistema de ecuaciones normales, no necesitamos aquí invertir la matriz cuadrada $A^TA$.
Procedemos a realizar la inversión. Partimos de $\mathbf{m}_0=(5,2)^T$, utilizamos $\epsilon=0.001$ y un número máximo de 100 iteraciones:
vector de parámetros inicial: a = 5.0 b = 2.0 vector de parámetros obtenido: a* = 10.051 (redondeo a tres decimales) b* = 4.931 (redondeo a tres decimales) número de iteraciones: 100 tamaño de paso: 0.001
Visualizamos el resultado del método gradient descent :
Eso es todo por hoy.
Esta actividad práctica está basada en la siguiente notebook; en el libro de deep learning y en esta gran referencia sobre métodos numéricos.