Regresión

Mínimos Cuadrados, Mìnimos Cuadrados No-Lineales y Logística

Mariano Rivera

febrero 2017

%matplotlib inline

Unas notas previas:
$\nabla_x \, y^\top x = y \;\;\;\;\;\;\;\;\;\; \\ \nabla_x \, A x = A \;\;\;\;\;\;\;\;\;\; \\ \nabla_x \, x^\top A x = A^\top x + A x \\ \;\;\;\;\;\;\;\;\;\;\;\; \nabla_x \, x^\top A x = 2 Ax, \;\; si \; A^\top=A \\ \;\;\;\;\;\;\;\;\;\;\;\; \nabla_x \, \|x-y\|^2 = \nabla_x \, (x-y)^\top (x-y) \\ \;\;\;\;\;\;\;\;\;\;\;\; \;\;\;\;\;\; = 2 x^\top (x-y) \\ \;\;\;\;\;\;\;\;\;\; \nabla^2_{xx} \, x^\top A x = 2A, \;\; si \; A^\top=A$

Regresión Lineal

Sea $\mathbf{x_i} \in \mathbb{R}^n$ un vector que corresponde al $i$-ésimo dato tal que la $j$-ésima entrada $\mathbf{x}_{ij}$ corresponde a la $j$-ésima característica. Luego arreglamos los datos de la forma

(1)
$X = \begin{bmatrix} \mathbf{x}^\top_1 \\ \mathbf{x}^\top_2 \\ \vdots \\ \mathbf{x}^\top_m \\ \end{bmatrix}, \\$
de tal que cada dato correponde a un renglón de la matrix $X \in \mathbb{R}^{m \times n}$.

Luego tenemos una observacion $\mathbf{y}_i \in \mathbb{R}^k$ (generalmente k=1) asociada a cada dato, variable dependiente. Y denotamos por

(2)
$Y = \begin{bmatrix} \mathbf{y}^\top_1 \\ \mathbf{y}^\top_2 \\ \vdots \\ \mathbf{y}^\top_m \\ \end{bmatrix},\\$
la matriz de variables dependientes, cadad renglón en una variable dependiente.

Regresión lineal multivariada

Ahora que asumimos que existencia de una función desconocida

(3)
$\mathbf{y}_i=f(\mathbf{x}_i) + \eta_i \\$
donde $f$ es “suave” y $\eta_i$ es ruido. El problema de regresión es estimar $f$ a partir de una serie de datos $X, Y$.
Existen muchas maneras de proponer la función $f$. Por ejemplo, usando polinomios:

(4)
$\mathbf{y}_{il} = \theta_{0l} + \theta_{1l} \mathbf{x}_{i1} + \theta_{2l}\mathbf{x}_{i2} + \ldots + \theta_{nl} \mathbf{x}_{in} + \eta_{li}$
La función (4) se pueded escribir como:

(5)
$Y = \mathbf{X} \Theta + \eta$

Donde hemos definido el vector

(6)
$\theta_l = \begin{bmatrix} \mathbf{\theta}_{0l} \\ \mathbf{\theta}_{1l} \\ \mathbf{\theta}_{2l} \\ \vdots \\ \mathbf{\theta}_{nl} \\ \end{bmatrix}, \;\;\; l=1,2,\ldots,k;$

Luego, definimos la matriz de parámetros

(7)
$\Theta = [\theta_1, \theta_2, \ldots , \theta_k ]$
donde la $l$-ésima columna corresponden a los parámetros para generar la componente $l$-ésima de $Y$. También hemos definido la matrix

(8)
$\mathbf{X} = [ \mathbf{1} \, | \, X]$
donde $X$ es definida en (1) y $\mathbf{1}$ es un vector con todas sus entradas igual a 1, cuyo tamaño depende del contexto.

Note que la función en (5) es lineal en los parámetros $\theta$

Es posible incorporar términos de mayor orden o no lineales en $\mathbf{x}$ y manter el modelo lineal respecto a los parámetros $\theta$:

$\mathbf{y}_{il} = \theta_{0l} + \theta_{1l} \mathbf{x}_{i1} + \theta_{2l}\mathbf{x}_{i2} + \ldots + \theta_{nl} \mathbf{x}_{in} + \\ \theta_{n+1,l} \mathbf{x}_{i,1}^2 + \theta_{n+2,l}\mathbf{x}_{i,2}^2 + \ldots + \theta_{2n,l} \mathbf{x}_{in}^2 + \\ \theta_{2n+1,l} \cos(\mathbf{x}_{i1}^2 +\mathbf{x}_{i2}^2 )+ \theta_{2n+2,l} \mathbf{x}_{i3} \mathbf{x}_{i4} + \theta_{2n+3,l} \ln (\mathbf{x}_{in}+1) + \eta_{li} \\$

En cuyo caso, lo que cambia es el cálculo de la matriz de diseño $\mathbf{X}$.

Mínimos Cuadrados Lineales

Una forma de resolver el problema de estimar los parámetros $\Theta$ de la función $f$ es mediante mínimos cuadrados lineales, capítulo 3 en [3]:
$\Theta^* = \underset{\Theta}{\operatorname{arg\,min}}\; F(\Theta) = \frac{1}{2} \| Y- \mathbf{X} \, \Theta \|_F^2 \\ (9)$
donde la norma de Frobenius esta definida por
$\| A\|_F^2 = \sum_i \sum_j a_{ij}^2.$
En este contexto $\mathbf{X}$ es llamada la “matriz de diseño”.

La función $F$ puede también ser escrita como

(10)
$F(\Theta) = \sum_l \sum_i \left[ Y_{il} - \sum_j \mathbf{X}_{ij} \Theta_{jl} \right]^2 = \sum_l \sum_i \left[Y_{il} - \mathbf{x}_i^\top\theta_{l} \right]^2 \\ \;\;\;\;\;\;\;\;\; =\sum_i \left[Y_{i1} - \mathbf{x}_i^\top\theta_{1} \right]^2 + \sum_i \left[Y_{i2} - \mathbf{x}_i^\top\theta_{2} \right]^2 + \ldots + \sum_i \left[Y_{ik} - \mathbf{x}_i^\top\theta_{k} \right]^2$
Es decir, el problema de estimación de mínimos cuadrados vector-valuado (la variable dependiente es un vector), se puede descomponer como múltiples problemas independientes de mínimos cuadrados univaluados (la variable dependiente es un escalar).

Solución Exacta al Problema de Mínimos Cuadrados Lineales

Como ya vimos, sin pérdida de generalidad, podemos concentrarnos solo en el caso univariado:

(11)
$\theta^* = \underset{\theta}{\arg \min}\; F(\Theta) = \frac{1}{2}\sum_i [Y_i - \mathbf{x}_i^T\theta]^2 = \frac{1}{2}\|Y- \mathbf{X}\, \theta \|_2^2$

Donde hemos incluido el factor $1/2$ por conveniencia en el algebra que desarrollaremos y dado que no afecta el problema:
$\underset{\theta}{\operatorname{arg\,min}}\; F(\theta) = \underset{\Theta}{\operatorname{arg\,min}}\; c F(\theta)$
para cualquier constante $c \in \mathbb{R}$.

Dado que (11) es convexa, el mínimo se encuetra resolviendo la condición de primer orden:

(12)
$\nabla_\theta F(\theta) = 0,$
esto es:

(13)
$\mathbf{X}^T\left(\mathbf{X} \theta - Y \right) =0$
Luego, $\theta$ se puede calcular resolviendo

(14)
$\mathbf{X}^\top\mathbf{X} \theta = \mathbf{X}^\top Y$
o directamente:

(15)
$\theta = \left( \mathbf{X}^\top\mathbf{X} \right)^{-1}\mathbf{X}^\top Y.$
Recordemos que en general $\mathbf{X}$ es singular dado que tenemos mas datos que parámetros a estimar $m > n+1$. Luego definimos la matriz
$M = \left( \mathbf{X}^\top\mathbf{X} \right)^{-1}\mathbf{X}^\top$
que es llamada la pseudo-inversa de Moore-Penrose de $\mathbf{X}$ dado que satisface $M \mathbf{X}=I$.

Si el problema es multivaluado, entonces:

(16)
$\theta_l = M Y_l; \;\; l=1,2, \ldots,k \\$
Es decir, todo se reduce a calular $M$.

Usando el modelo, podemos obtener la predicción de los valores de $Y_l$:

$\hat Y_l = \mathbf{X} \theta_l = P Y_l$
donde

$P = \mathbf{X} \left( \mathbf{X}^\top\mathbf{X} \right)^{-1}\mathbf{X}^\top$
es la “Matriz de Proyección” porque lleva los datos $Y_l$ al espacio generado (spanned) por el regresor. Lo que implica en si mismo un residual de modelado.

Solución iterativa al problema de mínimos cuadrados

Otra forma de resolver el problema de mínimos cuadrados es usando algorirmos iterativos, algo comun en problemas de optimización.

Esta estrategia de solución es preferida para problemas de gran escala: $m$'s o $n$'s grandes

Para ello usamos la formula de recurrencia

(17)
$\theta^{(t+1)} = \theta^{(t)} + \alpha^{(t)} p^{(t)}$
donde

1. $\theta^{(t)}$ es en valor de los parámetros en la iteración actual

2. $\theta^{(t+1)}$ es en valor de los parámetros actualizados

3. $p^{(t)}$ es una dirección de descenso. Garantiza que $F(\theta^{(t)} + \epsilon p^{(t)} ) \le F(\theta^{(t)})$ para una $\epsilon$ suficientemente pequeña.

4. $\alpha^{(t)}$ es el tamaño de paso

Nota: $p$ es una dirección de descenso ssi $p^\top \nabla F < 0$; la elección obvia es la denominada “dirección de máximo descenso”: $p= -\nabla F$, ¡Que no es necesariamente la mejor!

Obtengamos el gradiente parcial de F respecto a un parámetro, digamos el $\theta_ {\hat j}$:

(18)
$\frac{\partial F}{\partial \theta_{\hat j}} = \sum_i \left[Y_i - \mathbf{x}_i^T\theta \right] \frac{\partial }{\partial \theta_{\hat j}} \left[Y_i - \sum_j \mathbf{x}_{ij} \theta_j \right] \\ = \sum_i \left[Y_i - \mathbf{x}_i^T\theta \right] (-\mathbf{x}_{i{\hat j}} )$
Luego el gradiente esta dado por

(19)
$\nabla_\theta F(\theta) = \begin{bmatrix} \frac{\partial F}{\partial \theta_0} \\ \frac{\partial F}{\partial \theta_1} \\ \frac{\partial F}{\partial \theta_2} \\ \vdots \\ \frac{\partial F}{\partial \theta_n} \end{bmatrix} = - \begin{bmatrix} \sum_i \left[Y_i - \mathbf{x}_i^T \theta \right] 1\\ \sum_i \left[Y_i - \mathbf{x}_i^T \theta \right] \mathbf{x}_{i1} \\ \sum_i \left[Y_i - \mathbf{x}_i^T \theta \right] \mathbf{x}_{i2} \\ \vdots \\ \sum_i \left[Y_i - \mathbf{x}_i^T \theta \right] \mathbf{x}_{in} \end{bmatrix}$

Reducción del error en un sólo parámetro $\theta_{j}$

Podemos ver que en (19) se actualiza cada entrada de $\theta$ independientemente. Luego una entrada del vector de parámetros se actualiza con la regla

(20)
$\theta^{(t+1)}_j = \theta^{(t)}_j - \alpha \sum_i \left(\mathbf{x}_i^T \theta -Y_i \right) \mathbf{x}_{ij}$
donde hemos asumido que el tamaño de paso $\alpha$, conocido como “razón de aprendizaje” en el contexto de aprendizaje de máquina, es pequeño e igual para todas las iteraciones.

El procedimiento (20) es conocido como algoritmo de descenso de gradiente.

Reducción del error respecto a una sola muestra $\mathbf{x}_i$

Solo por curiosidad, tomamos la muestra $\mathbf{x}_i$ y ajustemos el $\theta_j$. Ahora la regla de aprendizaje queda

(21)
$\theta^{(t+1)}_j = \theta^{(t)}_j - \alpha \left( \mathbf{x}_i^\top \theta - Y_i \right) \mathbf{x}_{ij}$
que es es la regla de aprendizaje de Widrow-Hoff.

Se pueded observar que la magnitud de la actualización es proporcional al error $\left( Y_i - \mathbf{x}_i^\top \theta \right)$: si hay una muestra donde el error es pequeño, el ajuste será pequeño; si por otro lado para una muestra el error es grande, se hará un ajuste grande.

Ejemplo de Mínimos Cuadrados

Dado el modelo lineal
$y_i = \theta_0 + x_{i1} \theta_1 + x_2 \theta_{i2} + n_i$
Generaremos $m=10$ muestras y ajustamos los parámetros $\theta$ con:

1. Usando el solver en scipy.linalg.lsts

2. Pseudo-inversa de Moore-Penrose (numpy)

3. Descenso de gradiente

import numpy as np
import scipy.linalg as spla

m=100
n=2
theta = np.random.rand(3)
print(theta)

# Matriz de diseño
X = np.concatenate((np.ones((m,1)),np.random.rand(m,n)), axis=1)
n = 0.01*np.random.randn(m)

# variable dependiente
# y = np.dot(X,theta)+n
y = X@theta + n

[ 0.62450876  0.04722106  0.98380124]

Solver en scipy.linalg.lstsq

theta_ls=spla.lstsq(X,y)
theta_ls[0]

array([ 0.62408591,  0.04626653,  0.98832767])

Pseudoinversa de Moore-Penrose

# M = np.dot(la.inv(np.dot(X.T,X)),X.T)
# theta_mp=np.dot(M,y)

# con el nuevo operador de multiplicacion de matrices @ en vez de np.dot
theta_mp = (spla.inv(X.T@X)@X.T)@y
theta_mp

array([ 0.62408591,  0.04626653,  0.98832767])

Solución mediante descenso de gradiente simple

Formula de actualización iterativa:
$\theta_{t+1} = \theta_t + \alpha_t p_t$
donde

• $\theta_t$ es el valor actual de los parámetros

• $p_t$ es la dirección de descenso: $p_t \equiv -\nabla f(\theta_t)$ para el caso de descenso de gradiente.

• $\alpha_t$ es el tamaño de paso: $\alpha_t \equiv$ constante pequeña para decenso de gradiente simple.

def descensoLSTSQ(X, y, nIter=100, alpha=0.01):
    m,n = X.shape
    theta = np.random.rand(n)
    for t in range(nIter):
        theta = theta - alpha * (X@theta-y)@X  #np.dot((np.dot(X,theta)-y),X)       
    return theta

theta_dg=descensoLSTSQ(X,y, nIter=10000, alpha=0.001)
theta_dg

array([ 0.62408591,  0.04626653,  0.98832767])

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

plt.figure(figsize=(12,8))


fig = plt.figure(figsize=(12,8))
ax = fig.gca(projection='3d')

v = np.arange(0, 1, 0.1)
w = np.arange(0, 1, 0.1)
x1, x2 = np.meshgrid(v, w, sparse=True)

haty = theta_dg[0] + theta_dg[1]*x1+ theta_dg[2]*x2

ax.scatter(X[:, 1], X[:, 2], y, c='g')
ax.plot_surface(x1, x2, haty, rstride=8, cstride=8, alpha=0.4)
plt.show()

<Figure size 864x576 with 0 Axes>

Mínimos Cuadrados Lineales como un Problema de Estimación de Máxima Verosimilitud

Regresando al modelo mediante el cual asumimos se relacionan las variables dependientes $y$ y las independientes $x$:

(22)
$\mathbf{y}_i = \mathbf{x}_i^\top \theta + \eta_i.$
Donde habiamos dicho que $\eta_i$ era ruido. Bueno, este ruido tienen una distribución y puede estar, o no, correlacionado.

Asumamos que no es correlacionado y su distribución se mantienen no cambia dependiendo de las muestras; esto es es ruido Independiente e Identicamente Distribuido (IID).

Además asumamos, ahora, que tienen distribución Gaussiana con media cero y varianza $\sigma^2$: $\eta_i \sim \mathcal N (0, \sigma^2)$.

Luego, la densidad de $\eta_i$ esta dada por

(23)
$p(\eta_i) = \frac{1}{\sqrt{2\pi} \sigma} \exp \left[-\frac{\eta_i^2}{2\sigma^2} \right] \\$

Recordemos que $\int_{-\infty}^{\infty} p(\eta) d \eta =1$ pero la probabilidad de que $\eta_i=c$ esta dada por $\int_c^c p(\eta) \, d \eta =0$.

Luego de (22) obtenemos que $\eta_i = \mathbf{y}_i - \mathbf{x}_i^\top \theta$ lo que implica, usando (23), que

(24)
$p(\mathbf{y}_i | \mathbf{x}_i ; \theta) = \frac{1}{\sqrt{2\pi} \sigma} \exp \left[-\frac{1}{2\sigma^2} \left( \mathbf{y}_i - \mathbf{x}_i^\top \theta \right)^2 \right] \\$

$p(\mathbf{y}_i | \mathbf{x}_i ; \theta)$ se lee como la probabilidad condicional de $\mathbf{y}_i$ dado $\mathbf{x}_i$; usando los parámetros (no son variables aleatorias) $\theta$.

La notación $p(\mathbf{y}_i | \mathbf{x}_i ; \theta)$ indica que en nuestro modelo $ \mathbf{y}_i$ depende de $\mathbf{x}_i$ (variable dependiente e independiente, respectivamente) a través de los parámetros $\theta$.

Nosotros estamos interasados en estimar la función que relaciona los datos observados $\mathbf{y}$ y $\mathbf{X}$. Note que ya estamos hablando de todos dos datos, no solo del $i$-ésimo.

Esta función que relaciona a $\mathbf{y}$ y a $\mathbf{X}$ depende de los parámetros desconocidos $\theta$.

Por lo que es conveniente introducir una función que se vea explicitamente como función de $\theta$.

Función de verosimilitud

Esta función se denomina verosimilitud (likelihood):
(25)
$L(\theta) = L(\theta; \mathbf{y}, \mathbf{X}) = p( \mathbf{y} | \mathbf{X}; \theta)$

Donde $p( \mathbf{y} | \mathbf{X}; \theta)$ es la probabilidad de toda la muestra; (24) es la probabilidad del $i$-ésimo dato.

Dado que asumimos que $\eta_i$ es IDD:
(26)
$L(\theta) = p( \mathbf{y} | \mathbf{X}; \theta ) = \prod_i p(\mathbf{y}_i | \mathbf{x}_i ; \theta) \;\;\;\;\;\;\;\;\;\;\;\;\\ \;\;\;\;\;\;\;= \prod_i \frac{1}{\sqrt{2\pi} \sigma} \exp \left[-\frac{1}{2\sigma^2} \left( \mathbf{y}_i - \mathbf{x}_i^\top \theta \right)^2 \right] \\ \;\;\;\;\;\;\;= \frac{1}{\sqrt{2\pi} \sigma} \exp \left[-\frac{1}{2\sigma^2} \sum_i \left( \mathbf{y}_i - \mathbf{x}_i^\top \theta \right)^2 \right]$

Función de log-verosimilitud y los Mínimos Cuadrados

Dado que $\max_x f(x) = \max_x g(f(x))$ si $g$ es estrictamente creciente. Por ejemplo $g_1(x) = \log(x)$, $g_2(x) = x^2$ ambas cumplen para para $x>0$.

Luego en vez de maximizar la verosimilitud, uno pueded maximizar la log-verosimilitud: $\mathcal{l}(\theta) = \log L(\theta)$

(27)
$\mathcal{l}(\theta) = \log \frac{1}{(\sqrt{2\pi} \sigma)^n} \exp \left[-\frac{1}{2\sigma^2} \sum_i \left( \mathbf{y}_i - \mathbf{x}_i^\top \theta \right)^2 \right] \\ = \log \left[ \frac{1}{(\sqrt{2\pi} \sigma)^n} \right]-\frac{1}{2\sigma^2} \sum_i \left( \mathbf{y}_i - \mathbf{x}_i^\top \theta \right)^2$

Entonces, maximizar $\mathcal{l}(\theta)$ equivale a resolver la minimización:

(28)
$\min \frac{1}{2}\sum_i\left( \mathbf{y}_i - \mathbf{x}_i^\top \theta \right)^2$
que resulta ser el problema de mínimos cuadrados original.

¿Como llegamos a ello?

---

1. Asumimos un modelo generador que relacional las variables independientes a partir de las dependientes

2. Asumimos un modelo de ruido: IID y Gaussiana.

3. La suposición IID implicó suma de residuos (errores) independientes en cada muestra.

4. La distribución del ruido fue determinante que la log-verosimilitud se reduzca a la norma $L_2$

---

Nótese: otra distribución del ruido implicará otra norma.

Mínimos Cuadrados No-Lineales

Recordemos, los datos son denotados por

$\mathbf{x_i} \in \mathbb{R}^n \\ \mathbf{y_i} \in \mathbb{R} \\ X \in \mathbb{R}^{m \times n} \\ Y \in \mathbb{R}^{m}$

Como ya vimos, basta que analicemos la regresion monovaluada. Ahora, en el modelo (3) asumamos no-lineal función $f:\mathbb{R}^n \rightarrow \mathbb{R}$ y que dicha función es evaluada dados unos parámetros $\theta \in \mathbb{R}^k$:

(29)
$\mathbf{y}_i=f(\mathbf{x}_i; \theta) + \eta_i.\\$
Nóte que al escribir el modelo (29) hemos puesto un ‘;’ entre la $x$ y $\theta$. La razón es que queremos enfatizar que la función “produce” valores conforme cambiamos $x$, y que el hecho de que, por ahora, no conozcamos los valores de los parámetros $\theta$ es circunstancial. Una vez que la caractericemos (estimemos $\theta$) computaremos valores $f(x)$.

Para estimar $\theta$, definimos una función que depende explícitamente de su valor. Por ello, definimos el residual:

(30)
$r_i(\theta) = \mathbf{y}_i - f(\mathbf{x}_i, \theta)$

Luego, en un esquema de mínimos cuadrados (no lineales), los parámetros $\theta$ los puedemos estimar mediante la solución del problema de optmización:

(31)
$\theta^* = \underset{\theta}{\operatorname{arg\,min}}\; F(\theta) = \frac{1}{2}\sum_i r_i(\theta)^2 = \frac{1}{2} \| r(\theta) \|_2^2= \frac{1}{2} r(\theta)^\top r(\theta)\\$

donde hemos definido la función vectorial $r:\mathbb{R}^k \rightarrow \mathbb{R}^m$:

(32)
$r(\theta) = \begin{bmatrix} r_1(\theta) \\ r_2(\theta) \\ \vdots \\ r_m(\theta) \\ \end{bmatrix}$

Solución iterativa al probema de mínimos cuadrados no-lineales

Si usamos la estrategia de descenso de gradiente, entonces el gradiente esta dado por

(33)
$\nabla_\theta F(\theta) = \sum_i r_i(\theta) \nabla_\theta r_i(\theta)\;\;\;\;\;\;\;\;\;\;\;\; \\ \;\;\;\;\;\;\;\;\; = \nabla_{\theta} r(\theta) \; r(\theta) = J(\theta)^\top r(\theta)\\$

Donde definimos el “Jacobiano” $J$ de una función vectorial como:

(34)
$J_{ij} = \frac{\partial r_i }{\partial \theta_j}$
Entonces

(35)
$J(\theta)^\top = \nabla_\theta r(\theta) = \begin{bmatrix} \nabla_\theta r_1(\theta), \nabla_\theta r_2(\theta), \ldots, \nabla_\theta r_m(\theta) \end{bmatrix}$
Entonces, $J(\theta)^\top$ es el gradiente de cada una de las funciones $r_i$ acomodado columna por columna.

El algortimo de descenso de gradiente estará dado por:

(36)
$\theta^{t+1} = \theta^{t} - \alpha \sum_i r_i(\theta) \nabla_\theta r_i(\theta) \\ = \theta^{t} - \alpha \; J(\theta)^\top r(\theta) \;\;\;$

Newton, Gauss-Newton y Levenberg-Marquardt

Sin dar ninguna derivación formal, aqui introducimos directamente métodos mas eficientes de optimizacion de mínimos cuadrados no-lineales; ver su derivación en un texto de Optimizacion Numérica [2]_.

Los métodos de Newton calculan la actualización usando:

(37)
$\theta^{t+1} = \theta^{t} - \alpha \; p^{t}$

donde $\alpha$ se escoge tal que $F(\theta^{t+1}) \le F(\theta)$ y la dirección de descenso $p$ es solución de alguno de los tres siguientes:

1. Newton

(38)
$\left[ J(\theta)^\top J(\theta) + H(\theta) \, r(\theta) \right] p = - J(\theta)^\top r(\theta)$

donde
$H(\theta) r(\theta) = \sum_i r_i(\theta) \nabla_\theta^2 r_i(\theta)$
Tal que $H$ es el Hessiano de $r$ y $\nabla_\theta^2 r_i$ es el Hessiano de $r_i$.
2. Gauss-Newton

(39)
$J(\theta)^\top J(\theta) p = - J(\theta)^\top r(\theta)$

3. Levenberg-Marquardt (LM)

(40)
$\left[ J(\theta)^\top J(\theta) + \tau I \right] p = - J(\theta)^\top r(\theta)$
con $\tau > 0$ garantizando que la matriz tenga un número de condición apropiado para ser establemente resuelto el sistema.

El método de Newton no es el preferido debido a:

• Require calcular Hessianos costosos

• Si los residuales $r_i$ se distribuyen con media cero, en el óptimo es de esperarse que el término extra se cancele: suma de residos pesados

• Los datos atípicos (outliers), harían que la suma no se cancelaría e introducirían error en el método.

[2] J. Nocedal and S. J. Wright, Numerical Optimization, Springer 2nd Ed., 2006.

Gradiente estocástico

Como podemos observar la función de costo de mímimos cuadrados (lineales o no-lineales) puede escribirse como la suma de residuos individuales: la suma corre sobre todos los residuales. En el caso de problemas de aprendizaje de máquina, estos residuales se asocian con los datos y la suma se hace sobre toda la población. Esto admite una interpretación extra: la funcion de costo y su gradiente se pueden ver como esperanzas:

(41)
$F(\theta) = \frac{1}{2}\sum_{i=1}^m r_i(\theta)^2 = \frac{m}{2}\mathbb{E} \{ r_i(\theta)^2 \}$

(42)
$\nabla_\theta F(\theta) = \sum_{i=1}^m r_i(\theta) \nabla_\theta r_i(\theta) = m \,\mathbb{E} \{ r_i(\theta) \nabla_\theta r_i(\theta) \}$

Si en vez de considerar toda la población, consideramos a la media muestral como estimador de la función de costo y su gradiente:

(43)
$\hat F(\theta) = \frac{1}{2}\sum_{i \in \mathcal{S} } r_i(\theta)^2$

(44)
$\nabla_\theta \hat F(\theta) = \sum_{i \in \mathcal{S}} r_i(\theta) \nabla_\theta r_i(\theta)$
Donde $\mathcal{S} \subset \{1,2,3,\ldots,m \}$ determina la muestra. Si se usa el $\hat F(\theta)$ en vez de $F(\theta)$ en el procedimiento iterativo de optimización el algoritmo se denominará:

Gradiente estocástico, Newton-estocástico, Gauss-Newton estocástico, LM estocástico, Quasi-Newton estocástico, etc.

Estas variantes estocásticas de algoritmos deterministas han sido extensivamente usadas en problemas de gran escala y sobre todo en problemas relacionados con técnicas llamadas de APRENDIZAJE PROFUNDO (DEEP LEARNING).

Su ventaja es que evitan ser atrapados por mínimos locales y son computacionalmente eficientes.

Regresión Logística

Consideremos el problema de clasificación binaria donde los datos se denotan por

$\mathbf{x_i} \in \mathbb{R}^n \;\;\; \\ \mathbf{y_i} \in \{0, 1 \} \\ \;\; X \in \mathbb{R}^{m \times n} \;\;\; \\ \;\; Y \in \{0, 1 \}^{m}$

donde asumimos que a todos los datos les hemos agregado un $1$ en su primer entrada $x_{i1} = 1$ y a partir del 2do. elemento se tienen las características observadas. Esto simplifica la notación al escribir el sistema lineal.
Luego

$X = \begin{bmatrix} \mathbf{x}_1^\top \\ \mathbf{x}_2^\top \\ \mathbf{x}_3^\top \\ \vdots \\ \mathbf{x}_3^\top \end{bmatrix}$

El algoritmo de regresión logística es un algoritmo de clasificación binaria para clases linealmente separables, capítulo 4 en [3], pero donde existen datos que se confunden (no exactamente linelamente separables).

Es un poco extraño, es un clasificador, pero se le denomina “regresión” logística. Si somos un poco pacientes, entenderemos el porque del nombre.

Primero introducimos la función sigmoide (que tiene forma de sigma o ‘s’):

(45)
$\phi(z) = \frac{1}{1+exp(-z)}$

La función (45) se grafica en seguida usando Python

[3] T._Hastie et al., The elements of Statistical Learning, 2nd Ed. Springer, 2009.

%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt

def sigmoid(z):
    return (1/(1+np.exp(-z)))

T = 10
z = np.arange(-T,T, 0.1)

plt.figure(figsize=(10, 4))
plt.axvline(0.0,color='k')
plt.axhline(0.0,ls='dotted', color='k')
plt.axhline(1.0,ls='dotted', color='k')
plt.axhline(0.5,color='k')
plt.plot(z,sigmoid(z))
plt.ylim(-.1,1.1)
plt.xlabel('z')
plt.ylabel('$\phi(z)$')
plt.suptitle('Sigmoide')
plt.show()

Recordenos que el $i$-ésimo dato esta dado por el vector columna $\mathbf{x}_i$. Ahora, asumiremos que las etiquetas binarias estarán sujetas al modelo

(46)
$\mathbf{y}_i=\phi( \mathbf{x}_i^\top \omega ) + \eta_i$

donde los parámetros $\omega \in \mathbb{R}^k$ son tal que el hiperplano

$\mathbf{x}_i^\top \omega > 0 \;\; si \;\; \mathbf{y}_i=1$

y

$\mathbf{x}_i^\top \omega <0 \;\; si \;\; \mathbf{y}_i=0;$

luego, la sigmoide $\phi:\mathbb{R} \rightarrow [0,1]$ hará el resto del trabajo: ajustar los valores del hiperplano a las etiquetas. Los errores en classificación y residuales se representan por $\eta$.

Usando una estrategia de mínimos cuadrados no lineales, $\omega$ se puede estimar mediante:

(47)
$\arg\min _\omega F(\omega) = \frac{1}{2}\sum_{i=1}^m \left[\mathbf{y}_i-\phi(\mathbf{x}_i^\top \omega) \right]^2$

En este caso, el gradiente estará dado por las parciales

(48)
$\frac{ \partial F(\omega) }{\partial \omega_j} = \sum_{i=1}^m \left[\phi(\mathbf{x}_i^\top \omega) - \mathbf{y}_i \right] \phi(\mathbf{x}_i^\top \omega) \, [1-\phi(\mathbf{x}_i^\top \omega)] \, \mathbf{x}_{ij}$

donde usamos que $\frac{ \partial \phi(z) }{\partial z} = \phi(z)[1-\phi(z)]$.

Sin embargo, mejor tomemos un enfoque probabilista para derivar el algoritmo de ML para la Regresión Logística , y veamos si tenemos una ventaja.

Derivación Probabilística de la Regresión Logística

Asumamos $p$ la probabilidad de que ocurra un evento, luego $1-p$ será la probabilidad de que no ocurra.

Ahora, ¿que tanto es mas probable que ocurra el evento a que no ocurra? Para responder esta pregunta calculamos los momios (odd ratios), que es la razón:

(49)
$\frac{p}{1-p}$

que se interpreta como sigue. Si es valor toma un valor de 2, significa que es 2 veces mas probable que ocurra $p$ respecto a que no ocurra. Una forma mas cómoda de ver esta razón es tomado su logaritmo (natural), y definimos la función logit:

(50)
$\mathrm{logit}(p) = \log \frac{p}{1-p}$

Los momios y la función Logit son graficados usando Python con el código siguiente

%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt

def logit(p):
    return np.log(p/(1-p))

p = np.arange(0.01,1,0.01)
one = np.ones(p.shape)

plt.figure(figsize=(15, 8))
plt.suptitle('Momios')

plt.subplot(121)
plt.plot(p,p/(1-p))
plt.axvline(0.0, color='k')
plt.axvline(1.0,ls='dotted', color='k')
plt.axhline(0.0,color='k')
plt.xlabel('p')
plt.ylabel('razon')
#plt.grid(True)

plt.subplot(122)
plt.plot(p,logit(p))
plt.axvline(0.0, color='k')
plt.axvline(1.0,ls='dotted', color='k')
plt.axhline(0.0,color='k')
plt.xlabel('p')
plt.ylabel('logit')
#plt.grid(True)

plt.show()

Asi que si $p=2/3$, tendiamos que $1-p = 1/3$ los momios son $2$.

Si $p=.9$, tendiamos que $1-p = .1$ los momios son $9$.

Si $p=.95$, tendiamos que $1-p = .05$ los momios son $19$.

Por ello, la gráfica anterior de la izquierda se dispara para valores de p cercanos a $1$. Por otro lado $\mathrm{logit}(p) = \log p - \log(1-p)$ es más fácil de entender, es cero para cuando $p=1-p=.5$, y antisimérica a partir de ese punto.

Verosimilitud de la Regresión Logística

Denotamos por $p(y_i=1|\mathbf{x}_i;\omega)$ la probabilidad condicional de clasificar el dato $i$ a la clase 1 dado el vector de rasgos $\mathbf{x}_i$, asumidos los parámetros $\omega$.

Luego, para cualquier clase:

(51)
$p(y_i|\mathbf{x}_i;\omega) = \left\{ \begin{matrix} \phi(\mathbf{x}_i^\top \omega) & y_i=1 \\ 1-\phi(\mathbf{x}_i^\top \omega) & y_i=0 \\ \end{matrix}\right.$

o podemos usar

(52)
$p(y_i|\mathbf{x}_i;\omega) = y_i \phi(\mathbf{x}_i^\top \omega) + (1-y_i)[1-\phi(\mathbf{x}_i^\top \omega)]$
o, también

(53)
$p(y_i|\mathbf{x}_i;\omega) = \phi(\mathbf{x}_i^\top \omega)^{y_i} [1-\phi(\mathbf{x}_i^\top \omega)]^{(1-y_i)}$
Las tres son equivalentes. Sin embargo, usaremos la forma (53) que corresponde a la distribución Bernoulli. Una distribución dicotómica, con solo dos posibles resultados: la probabilidad de éxitos $p$ y fracasos $(1-p)$.

Ahora asumamos la probabibilidad condicional de las clasificación de toda la muestra

(54)
$P(\mathbf{Y}|\mathbf{X};\omega) = \prod_i p(y_i|\mathbf{x}_i;\omega) = \prod_i \phi(\mathbf{x}_i^\top \omega)^{y_i} [1-\phi(\mathbf{x}_i^\top \omega)]^{(1-y_i)}$

donde asumimos que los datos estan IID.

Ahora, notamos que la condicional $P(\mathbf{Y}|\mathbf{X};\omega)$ es función de $\mathbf{Y}$, asume que ensayamos $\mathbf{X}$ para unos parámetros $\omega$. Pero si lo que queremos es aprender es $\omega$ a partir de datos de entrenamiento $[\mathbf{X}, \mathbf{Y}]$, necesitamos definir una función que depende explícitamente de $\omega$ y use como datos a $[\mathbf{X}, \mathbf{Y}]$. Esta función es, como antes, la verosimilitud:

(55)
$L(\omega; \mathbf{X}, \mathbf{Y} ) \overset{def}{=} P(\mathbf{Y}|\mathbf{X};\omega)$

o simplemente $L(\omega)$.

La log-verosimulitud

Para estimar $\omega$ debe podemos usar (55) como función de mérito. Es decir, encontrar la $\omega$ que maximice la verosimilitud. Pero esta función tiene productos de probabilidades (valores entre cero y uno) y podemos tener el riesgo de sobreflujo en los cálculos. Asi que mejor calculamos la log-verosilitud :
$l(\omega) = \log L(\omega)$
y tenemos:

(56)
$l(\omega) = \sum_i y_i \log \phi(\mathbf{x}_i^\top \omega) + (1-y_i) \log [1-\phi(\mathbf{x}_i^\top \omega)]$

Luego, para calcular el gradiente, calculamos cada parcial como

(57)
$\frac{\partial l(\omega)}{\partial \omega_j} = \sum_i \frac{\partial}{\partial z_i} \left\{ y_i \log \phi(z_i) + (1-y_i) \log [1-\phi(z_i)] \right\} \frac{\partial \phi(z_i) }{\partial \omega_j} \\ = \sum_i \left[ y_i - \phi(\mathbf{x}_i^\top \omega) \right] \mathbf{x}_{ij} \;\;\;\;\;\; \;\;\;\;\;\; \;\;\;\;\;\; \;\;\;\;\;\;$

Donde hemos usado

(58)
$\frac{\partial \phi(z)}{\partial z} = \frac{\partial}{\partial z} \frac{1}{1+ e^{-z}} = \frac{e^{-z}}{(1+ e^{-z})^2} \\ = \left(\frac{1}{1+ e^{-z}}\right)\left( \frac{1-1 + e^{-z}}{1+ e^{-z}}\right) \\ = \phi(z)[1- \phi(z)] \;\;\;\;\;\;$

(59)
$\frac{\partial}{\partial z} \left\{ y \log \phi(z) + (1-y) \log [1-\phi(z)] \right\} = \left(\frac{y}{\phi(z)} - \frac{1-y}{1-\phi(z)} \right) \frac{\partial \phi(z)}{\partial z} \\ = \left(\frac{y}{\phi(z)} - \frac{1-y}{1-\phi(z)} \right) \phi(z){(1- \phi(z)} \\ = y [1-\phi(z)] -(1-y) \phi(z) \\ = y - \phi(z)$

El gradiente de la función de de costo logística (56) esta dada por (57), que es muy similar en forma al gradiente de la regresión lineal de mínimos cuadrados, ver (18). Por lo computacionalmente no es mucho mas costosa la regresión logística.

La regresión logistica puede explicarse a partir de la siguiente figura.

1. Se tienen dos clase (casi) linealmente separables (puntos rojos y azules), con etiquetas binaria $y_i \in \{-1, 1\}$.

2. El plano $\omega^\top x$ es tal que al evaluarse en cualquier punto $x_i$ con etiqueta $y_i =-1$ resulta en $\omega^\top x_i <0$. Similarmente, para $x_j$ con etiqueta $y_j=-1$ resulta en $\omega^\top x_j >0$. Es importante notar que existen algunos puntos que no cumplen con la regla que hemos establecido (en la figura hay un punto rojo), el propósito es calcular $\omega$ tak que minimize los errores.

3. Al aplicar la función sigmoide al plano, se tienen que los puntos el mapeo queda en el intervalo $[-1, 1]$, y se aproximará a la saturación ($\{-1, 1\}$) a medida que la pendiente sea máxima: en el límite, es el plano vertical.

Es decir, la optimización de la regresión logística no esta acotada si las clases son perfectamente separables. Ésto se resuelve combiando al azar la etiquerta de un dato (convirtiendolo en outlier).

¿Que pasa si en vez de construir nuestro regresor sobre la restricción $y_i(\omega^\top x_i) >0$, construimos sobre $y_i(\omega^\top x_i) > 1$?

Respuesta. Ver las Máquinas de Vectores de Soporte [4]

[4] C. Cortes and V. Vapnik. “Support-vector networks”. Machine Learning. 20 (3): 273–297 (1995).