(769) - Expresiones para soluciones de ecuaciones diferenciales. ¿Existen o se crean?

Poco antes de verano estaba con un amigo del doble grado y redactor de este blog hablando de ecuaciones diferenciales y salieron dos ejemplos en particular: la del péndulo realista, $y^{\prime\prime}+\beta\,y^\prime+{\omega_0}^2\sin(y)=0$ (con amortiguamiento/fricción y sin aproximaciones de ángulo pequeño), y el sistema de Lotka-Volterra: $\displaystyle \begin{cases} x^\prime & = & x(\alpha-\beta y) \\ y^\prime & = & -y(\gamma-\delta x) \end{cases}$ . Él me preguntó que si acaso tenían solución, y yo le dije que sí, que se podían construir. Él realmente no me inquirió en si eran resolubles (ya fueren analítica o numéricamente), sino si existía alguna función conocida, que pudiese ser composición y/o combinación de otras funciones elementales o no (por ejemplo) tal que fuese solución suya, es decir, si se podía expresar su solución en términos de funciones conocidas. Él se quejó de mi respuesta, a lo que yo le dije que realmente eso poco importaba a la hora de resolver ecuaciones diferenciales y al estudiarlas. A veces la propia solución puede ser tan engorrosa analíticamente que es reconfortante pensar que siempre nos quedará el cálculo numérico...

Por ejemplo, tomemos la función error gaussiano, una función no-elemental, que se define como: $\displaystyle \operatorname{erf}(x)\overset{\text{def}}{=}\frac{2}{\sqrt{\pi\;}}\int_0^x e^{-t^2}\;\text{d}t$ . Resolvamos la ecuación diferencial ordinaria de I orden $y^\prime = e^{-x^2}$ . Ahora bien, ¿la solución de esta ecuación diferencial es debida a dicha definición anterior, tal que $\displaystyle y\triangleq\frac{\sqrt{\pi\;}}{2}\operatorname{erf}(x)+y_0$ , o bien la definición anterior se dio para que satisficiere la ecuación diferencial $y^\prime = e^{-x^2}$ ?, es decir, ¿qué vino antes: la definición (y después la ecuación diferencial que resulta ser que la tiene como solución salvo constantes) o la ecuación diferencial (y la definición se dio posteriormente como su solución a este problema)?

Sin la definición de $\operatorname{erf}(x)$, el problema es irresoluble analíticamente, pero con dicha definición es trivial, ya que la propia definición es la solución. Es como si hubiésemos definido que esta función sea su solución. Es equivalente a construir la solución, y darla como definición, por lo que el problema ya está resuelto.

Esto puede parecer una tontería, o un poco enrevesado, pero si se quiere resolver $y^\prime = 3x^2$ , que tiene como solución $y=x^3+y_0$ , no se requiere de ninguna definición previa que no esté contemplada ya.

Uno siempre puede decir: Sea $y\overset{\text{def}}{=}\theta(t)=\theta(t|\beta,\omega_0;y_0,{y_0}^\prime)$ la solución a la ecuación diferencial $y^{\prime\prime}+\beta\,y^\prime+{\omega_0}^2\sin(y)=0$ en la variable $t$ donde $\beta,\omega_0\geqslant 0$ son parámetros que dependen de la ecuación diferencial, con $y_0,{y_0}^\prime$ como condiciones iniciales (hay soluciones analíticas si se linealiza o si al menos uno de los parámetros se anula, aunque con $\beta=0$ es enrevesada). Después mediante cálculo y análisis numéricos, con un conocimiento previo de ecuaciones diferenciales, uno puede estudiar la solución numéricamente, y a veces a lo sumo, hasta cierto punto analíticamente.

Otro ejemplo puede ser el oscilador de Van der Pol: $y^{\prime\prime}-\mu(1-y^2)y^\prime+{\omega_0}^2y=0$ (nótese que el caso $\mu=0$ es el oscilador armónico simple, por lo que sí tiene una solución analítica en este caso), otra vez ahí definir $y\overset{\text{def}}{=}\operatorname{vdp}(t)=\operatorname{vdp}(t|\mu,\omega_0;y_0,{y_0}^\prime)$ como una solución obtenida numéricamente y a partir de ella estudiar la solución "analítica".

Lo que estoy argumentando es equivalente a si un matemático en el $\text{siglo XVI-XVII}$, desconociendo tanto el número $e$ como las funciones exponencial $e^x$ y logaritmo $\ln(x)$ , decidiese analizar la ecuación diferencial ordinaria de primer orden $y^\prime=k\, y$ definiendo su solución como $y\overset{\text{def}}{=}\operatorname{sol}(t) = \operatorname{sol}(t|k;y_0)$ y estudiarla a partir de su solución construida. Hoy sabemos que tiene una forma cerrada, $y=y_0\, e^{k(t-t_0)}$ .

Uno puede incluso argumentar que estos ejemplos están sesgados y que no tienen validez, pero no es así, hay muchas ecuaciones diferenciales que no tienen necesariamente una solución analítica y para entenderlas es necesario su estudio con cálculo numérico y definir sus soluciones. En especial, en física, como dijo Steven Strogatz ($1959-$): "Desde Newton, la humanidad se ha dado cuenta que las leyes de la física siempre se expresan en la lengua de ecuaciones diferenciales." Veamos algunos ejemplos:

El sistema de $2$ ecuaciones diferenciales ordinarias de II orden que rigen el péndulo doble $$ \begin{cases} \ddot{\theta}_1 & = & \displaystyle \frac{-g (2m_1+m_2)\sin(\theta_1)-m_2g\sin(\theta_1-2\theta_2)-2\sin(\theta_1-\theta_2)m_2\big({\dot{\theta_2}}^2 l_2+{\dot{\theta_1}}^2 l_1\cos(\theta_1-\theta_2)\big)} {l_1\big(2m_1+m_2-m_2\cos(2\theta_1-2\theta_2)\big)} \\ \ddot{\theta}_2 & = & \displaystyle \frac{2 \sin(\theta_1-\theta_2) \big({\dot{\theta}_1}^2 l_1 (m_1 + m_2)+ g(m_1 + m_2) \cos(\theta_1) + {\dot{\theta}_2}^2 l_2 m_2 \cos(\theta_1 - \theta_2)\big)} {l_2 \big(2 m_1 + m_2 - m_2 \cos(2 \theta_1 - 2 \theta_2)\big)} \end{cases} $$ El sistema de $2$ ecuaciones diferenciales ordinarias de II orden que rigen el péndulo elástico $$ \begin{cases} \ddot{r} & = & \displaystyle (r+l_0)\dot{\theta}^2 -\frac{k}{m}r + g\cos(\theta) \\ \ddot{\theta} & = & \displaystyle -\frac{g\sin(\theta)+2\dot{r}\dot{\theta}}{r+l_0} \end{cases} $$ La ecuación diferencial ordinaria de II orden que rige el péndulo esférico $$ \ddot{\theta}(t)- \left(\frac{p_\varphi}{ml^2}\right)^2\frac{\cos\big(\theta(t)\big)}{\sin^3\big(\theta(t)\big)}+\frac{g}{l}\sin\big(\theta(t)\big) = 0 $$ El sistema de $6$ ecuaciones diferenciales ordinarias de II orden que rigen el problema de los dos cuerpos $$ \begin{cases} \ddot{\vec{r}}_1 & = &\displaystyle Gm_2 \frac{\vec{r}_2-\vec{r}_1}{{\|\vec{r}_2-\vec{r}_1\|_2}^3}\\ \ddot{\vec{r}}_2 & = &\displaystyle Gm_1 \frac{\vec{r}_1-\vec{r}_2}{{\|\vec{r}_1-\vec{r}_2\|_2}^3} \end{cases} \implies m_1\ddot{\vec{r}}_1 + m_2\ddot{\vec{r}}_2 = \vec{0} $$ El sistema de $3$ ecuaciones diferenciales ordinarias de I orden conocido como el atractor de Lorentz, un ejemplo famoso de la Teoría del Caos: $$ \begin{cases} \dot{x} & = & \sigma (y-x) \\ \dot{y} & = & x (\rho-z)-y \\ \dot{z} & = & xy-\beta z \end{cases} $$
Autor: Đɑvɪẟ Ƒernández-De la Cruʒ.

(761) - Mínimos cuadrados en machine learning. Least trimmed squares (LTS)

En el día de hoy veamos cómo se utiliza la regresión lineal en machine learning. Recapitulemos un poco cómo funcionan los mínimos cuadrados ordinarios, ordinary least squares u OLS por sus siglas en inglés, el caso más simple: Dadas dos sucesiones de valores, una predictora $\{x_i\}_{i=1}^n$ (abscisas), y otra respuesta $\{y_i\}_{i=1}^n$ (ordenadas), se buscan los argumentos (parámetros) que minimicen la siguiente restricción $Z_{OLS}$ (del alemán Zwang - obligación): $$ Z_{OLS} = \sum_{i=1}^n \big(y_i-(\beta_0+\beta_1x_i)\big)^2 \qquad \hat{\beta}_0,\hat{\beta}_1 = \operatorname{arg\,min} Z_{OLS}$$ Nótese que esta restricción siempre es no-negativa, $Z_{OLS}\geqslant 0$ , y la igualdad (que sería lo deseable) solo se da si el ajuste es perfecto y no hay ninguna desviación.

Dado que queremos hallar la recta que mejor se ajuste a $\{y_i\}_{i=1}^n$ en la expresión aparece $\beta_0+\beta_1x_i$ , sino, si supiésemos o si intuyésemos que los $\{y_i\}_{i=1}^n$ siguen una función $f(x)$ genérica (polinómica, sinusoidal, exponencial,...), entonces se buscaría qué parámetros de $f(x)$ minimizan $\displaystyle Z_{OLS} = \sum_{i=1}^n \big(y_i-f(x_i)\big)^2$ . Se suele denotar por $\hat{y}_i=f(x_i)$ a los valores predichos. El método de mínimos cuadrados busca los argumentos, i.e. parámetros, $\hat{\beta}_0,\hat{\beta}_1,\cdots$ que minimizan la suma total de los cuadrados de todas las $n$ diferencias entre los valores reales $y_i$ , y los valores predichos $\hat{y}_i$ , $Z_{OLS}$ , (esta diferencia es el residual $r_i\overset{\text{def}}{=}y_i-\hat{y}_i$, por eso realmente se suele poner como hallar los argumentos que minimizan $\displaystyle Z_{OLS}=\sum_{i=1}^n {r_i}^2$ ). En nuestro caso para la regresión de una recta si denotamos como $\hat{\beta}_0,\hat{\beta}_1$ a la ordenada en el origen y pendiente predichas según el modelo respetivamente, se tiene que $\hat{y}_i = \hat{\beta}_0+\hat{\beta}_1x_i$ .

Sin embargo, en machine learning no se utiliza esto, al menos no así. Se hace lo que se conoce como mínimos cuadrados recortados, least trimmed squares o LTS por sus siglas en inglés: Primero se preselecciona un número $k$ de puntos que se tendrán en cuenta (y los $(n-k)$ restantes para un estudio posterior), donde $k$ suele estar entre el $75\%-80\%$ del total $n$ de puntos. Empero dichos $k$ puntos que se toman no son al azar, sino que se escogen de una manera muy particular: de entre todos los $n$ puntos $\big\{(x_i,y_i)\big\}_{i=1}^n$ , se busca entre todas las $\displaystyle \binom{n}{k}$ posibles combinaciones (de $k$ elementos tomados de un supraconjunto de $n$ elementos) aquella cuya restricción $Z$ de $k$ sumandos sea la mínima de todas las posibles, es decir, en cualquier otro subconjunto de tamaño $k$ sobre el supraconjunto $n$ se ajustan peor entre sí a mínimos cuadrados (su restricción $Z$ es mayor).

Sea $\mathcal{Z}$ el conjunto de las $\displaystyle |\mathcal{Z}|=\binom{n}{k}$ posibles mínimas sumas parciales de los cuadrados de $k$ residuales, $\displaystyle \mathcal{Z} \overset{\text{def}}{=} \bigg\{ \min\!\sum_{i=1}^k {r_{\sigma(i)}}^2 \;|\; k\leqslant n\, ,\,\sigma\in S_n\bigg\}$ con $k$ fijo , donde $\sigma(\cdot)$ indica una permutación de los $n$ puntos. Buscamos qué argumentos minimizan $Z_{LTS}\overset{\text{def}}{=}\min\{\mathcal{Z}\}\geqslant 0$ . Por esta definición el coeficiente de correlación $R^2$ de $Z_{LTS}$ es el máximo de los del conjunto $\mathcal{Z}$ , es decir, $R^2(Z_{LTS})\triangleq\max\big\{R^2(\mathcal{Z})\big\}$ .

Es más, si denotamos por $\big\{{r_{(i)}}^2\big\}_{i=1}^n$ la sucesión de los cuadrados de los residuales ordenados de menor a mayor, es decir, ${r_{(i)}}^2\leqslant {r_{(i+1)}}^2$ , buscamos qué parámetros minimizan la suma parcial $\displaystyle Z_{LTS}=\sum_{i=1}^k {r_{(i)}}^2$ , que es menor o igual que la suma total $\displaystyle \sum_{i=1}^n {r_{(i)}}^2=\sum_{i=1}^n {r_i}^2 = Z_{OLS}$ . Esta última restricción es como ya hemos dicho, la que se busca minimizar en OLS. La suma parcial será igual a la suma total si y solo si $r_{(i)}=0 \quad i = 1,\cdots, n$, es decir, si todos, los $n$ puntos, están perfectamente alineados.

Para comparar ambos resultados, ya hemos visto que $0\leqslant Z_{LTS}\leqslant Z_{OLS}$ , pero estaría bien comparar sendos errores cuadráticos medios, es decir, $\displaystyle \frac{1}{k} Z_{LTS}$ frente a $\displaystyle \frac{1}{n} Z_{OLS}$ , o también $\displaystyle \frac{1}{k-p} Z_{LTS}$ frente a $\displaystyle \frac{1}{n-p} Z_{OLS}$ (que se pueden entender como estimadores de las varianzas de sendos residuales) donde $p$ es el número de parámetros que se estiman (en una recta por ejemplo son dos: la ordenada en el origen y la pendiente), siendo $(k-p)$ y $(n-p)$ sendos grados de libertad.

A este conjunto de $k$ elementos que minimizan la suma parcial de los $k$ cuadrados de los residuales entre todas las posibles se lo conoce como conjunto de entrenamiento (training data set en inglés), mientras que el conjunto restante de $(n-k)$ elementos se lo conoce como conjunto de validación (test data set en inglés). El porqué de separar estos dos es debido a que el conjunto con los $(n-k)$ elementos restantes puede estar corrupto, contaminado, o ser defectuoso (son los outliers). Por eso este método es robusto: considera el conjunto de $k$ elementos que mejor se ajustan entre sí, mientras que el conjunto de $(n-k)$ elementos se utiliza para comprobar cómo de bueno es nuestro modelo, si las hipótesis son las correctas, y si hay que replantearse algo. Sin embargo, el conjunto de $(n-k)$ elementos no influyen en hallar los argumentos que minimizan la restricción $Z_{LTS}$ .

Cabe resaltar que $\displaystyle \frac{n^k}{k!}\geqslant \binom{n}{k} \geqslant \frac{(n-k+1)^k}{k!}$ siendo equivalentes cuando $n\to\infty$, es decir, que aumentando $n$ (al aumentar también $k$ ) implica que las posibles elecciones tienen un crecimiento exponencial-potencial, por lo que hallar qué conjunto de $k$ elementos es el adecuado puede ser muy costoso computacionalmente, en especial con $n$ medianos o con algoritmos no-óptimos. Por otro lado intentar estimar el conjunto de los $k$ con una simulación tipo Monte Carlo no parece ser muy fructífera, ya que la probabilidad de acertar aleatoriamente es $\displaystyle \frac{1}{\binom{n}{k}}$ .

Como último comentario, también existe un método que generaliza OLS conocido como mínimos cuadrados ponderados, weighted least squares o WLS por sus siglas en inglés, donde en vez de minimizar $\displaystyle Z_{OLS}=\sum_{i=1}^n {r_i}^2$ , se busca minimizar $\displaystyle Z_{WLS}=\sum_{i=1}^n \omega_i\,{r_i}^2$ donde $\{\omega_i\}_{i=1}^n$ es una sucesión de pesos no-negativos con cada $\omega_i$ asociado a su respectivo $r_i$ (realmente a $x_i$ , pero es equivalente) y normalmente normalizados, $\displaystyle \sum_{i=1}^n \omega_i=1$ . Nuestro modelo de LTS es en realidad un caso muy particular de WLS donde $\displaystyle\omega_{\sigma(i)}=\begin{cases} \displaystyle \frac{1}{k} & i\in\{1,\cdots,k\} \\ 0 & i\in\{k+1,\cdots,n\} \end{cases}$ . Es más, se podría generalizar esto aún más a un modelo de mínimos cuadrados ponderados y recortados, weighted least trimmed squares o WLTS por sus siglas en inglés donde se buscan: $$\operatorname{arg\,\min} \sum_{i=1}^k \omega_{(i)}\,{r_{(i)}}^2 $$
Autor: Đɑvɪẟ Ƒernández-De la Cruʒ.