Notas para contenedor de docker:
Comando de docker para ejecución de la nota de forma local:
nota: cambiar <ruta a mi directorio> por la ruta de directorio que se desea mapear a /datos dentro del contenedor de docker.
docker run --rm -v <ruta a mi directorio>:/datos --name jupyterlab_r_kernel_tidyverse -p 8888:8888 -d palmoreck/jupyterlab_r_kernel_tidyverse:1.1.0password para jupyterlab: qwerty
Detener el contenedor de docker:
docker stop jupyterlab_r_kernel_tidyverseDocumentación de la imagen de docker palmoreck/jupyterlab_r_kernel_tidyverse:1.1.0 en liga.
Nota generada a partir de liga
En esta nota revisamos a los mínimos cuadrados lineales con y sin regularización. La regularización que utilizamos es la de Tikhonov también nombrada $\ell_2$ o ridge y la $\ell_1$ o también conocida como lasso (least absolute shrinkage and selection operator, Tibshirani, 1996). Se muestra el uso de métodos de descenso (ver 4.2.Algoritmos_para_UCO) para resolver los problemas de optimización que surgen en los modelos anteriores y no se tiene por objetivo la interpretación de los coeficientes estimados. Se comparan los resultados del paquete glmnet stanford, glmnet cran de R con los obtenidos en la implementación hecha por el prof en algoritmos/R, en específico algoritmos/R/algorithms_for_uco.R para problemas tipo UCO (Unconstrained Convex Optimization).
Obsérvese que hay una gran cantidad de modelos por mínimos cuadrados, por ejemplo:
Lineales u ordinarios (nombre más usado en Estadística y Econometría).
Cada uno de los modelos anteriores tienen diversas aplicaciones y propósitos. Los lineales son un caso particular del problema más general de aproximación por normas:
$$\displaystyle \min_{x \in \mathbb{R}^n} ||Ax-b||$$donde: $A \in \mathbb{R}^{m \times n}$, $b \in \mathbb{R}^m$ son datos del problema, $x \in \mathbb{R}^n$ es la variable de optimización y $|| \cdot||$ es una norma en $\mathbb{R}^m$.
Se asume en toda la nota que $m \geq n $ (más renglones que columnas en $A$).
$x^* = \text{argmin}_{x \in \mathbb{R}^n} ||Ax-b||$ se le nombra solución aproximada de $Ax \approx b$ en la norma $|| \cdot ||$. El vector: $r(x) = Ax -b$ se le nombra residual del problema.
Comentario: el problema de aproximación por normas también se le nombra problema de regresión. En este contexto, los vectores $a_1, a_2, \dots, a_n$ (columnas de $A$) son nombradas regresoras o vector de atributos y el vector $\displaystyle \sum_{j=1}^n x_j^*a_j$ con $x^*$ óptimo del problema es nombrado la regresión de $b$ sobre las regresoras. $b$ es la respuesta.
Si en el problema de aproximación de normas anterior se utiliza la norma Euclidiana o norma $2$, $|| \cdot ||_2$, y se eleva al cuadrado la función objetivo se tiene:
$$\displaystyle \min_{x \in \mathbb{R}^n} ||Ax-b||^2_2$$que es el modelo por mínimos cuadrados lineales cuyo objetivo es minimizar la suma de cuadrados de las componentes del residual $r(x)$.
A partir de aquí, la variable de optimización será $\beta$ y no $x$:
Supóngase que se han realizado mediciones de un fenómeno de interés en diferentes puntos $x_i$'s resultando en cantidades $y_i$'s $\forall i=0,1,\dots, m$ (se tienen $m+1$ puntos) y además las $y_i$'s contienen un ruido aleatorio causado por errores de medición:
El objetivo de los mínimos cuadrados es construir una curva, $f(x|\beta)$ que "mejor" se ajuste a los datos $(x_i,y_i)$, $\forall i=0,1,\dots,m$. El término de "mejor" se refiere a que la suma: $$\displaystyle \sum_{i=0}^m (y_i -f(x_i|\beta))^2$$ sea lo más pequeña posible, esto es, a que la suma de las distancias verticales entre $y_i$ y $f(x_i|\beta)$ $\forall i=0,1,\dots,m$ al cuadrado sea mínima. Por ejemplo:
Obs:
Si $m=3$ y $A \in \mathbb{R}^{3 \times 2}$ geométricamente el problema de mínimos cuadrados lineales se puede visualizar con el siguiente dibujo:
donde: $r(\beta) = y-A\beta$, el vector $y \in \mathbb{R}^m$ contiene las entradas $y_i$'s y la matriz $A \in \mathbb{R}^{m \times n}$ contiene a las entradas $x_i$'s o funciones de éstas $\forall i=0,1,\dots,m$.. Por el dibujo se tiene que cumplir que $A^Tr(\beta)=0$, esto es: las columnas de $A$ son ortogonales a $r(\beta)$. La condición anterior conduce a las ecuaciones normales:
$$0=A^Tr(\beta)=A^T(y-A\beta)=A^Ty-A^TA\beta.$$Finalmente, considerando la variable de optimización $\beta$ y al vector $y$ tenemos: $A^TA \beta = A^Ty$.
con $A[i,:]$ $i$-ésimo renglón de $A$ visto como un vector en $\mathbb{R}^n$. Es común dividir por $2$ la función objetivo para finalmente tener el problema:
$$\displaystyle \min_{\beta \in \mathbb{R}^n} \quad \frac{1}{2}y^Ty-\beta^TA^Ty + \frac{1}{2}\beta^TA^TA\beta.$$En cualquier reescritura de la función $f_o$, el problema de aproximación con normas, o bien en su caso particular de mínimos cuadrados, es un problema de optimización convexa (ver 4.1.Optimizacion_numerica_y_machine_learning).
En lo que sigue se utiliza una nomenclatura similar del paquete glmnet de R.
Una técnica muy utilizada en el contexto de machine learning es la regularización, la cual tiene diferentes efectos en la solución de los problemas que surgen en esta área (por ejemplo lidiar con multicolinealidad entre variables, ver Multicollinearity, o el sobre ajuste, ver Overfitting) . La regularización es un caso particular del problema más general de optimización multicriterio, multiobjetivo, vectorial o también nombrada Pareto, ver Multi objective optimization.
Al añadir regularización al problema de aproximación por normas, se obtiene un problema de optimización bi criterion en el que además de minimizar la norma $||A\beta-y||$, se tiene que encontrar $\beta \in \mathbb{R}^n$ con norma $||\cdot||$ lo más pequeña posible. Esto es, se debe resolver el siguiente problema de optimización convexa con dos objetivos $||A\beta-y||$, $||\beta||$:
$$\displaystyle \min (||A\beta-y||,||\beta||)$$respecto a $\mathbb{R}^2_+ = \{(u,v) \in \mathbb{R}^2 : u \geq 0, v \geq 0\}$.
Comentario: en este problema se tiene el tradeoff entre tener $||\beta||$ mínima y $||A\beta-y||$ "grande" o mínima $||A\beta-y||$ y $||\beta||$ "grande".
La regularización es una técnica para resolver el problema anterior pues se propone una función objetivo como una suma ponderada de los dos objetivos anteriores:
$$\displaystyle \min_{\beta \in \mathbb{R}^n} ||A\beta-y|| + \lambda ||\beta||$$donde: $\lambda > 0 $ es un parámetro del problema. En esta formulación $\lambda$ varía en $(0, \infty)$ y permite realizar el tradeoff en el tamaño entre $||A\beta-y||$ vs $||\beta||$ descrito anteriormente.
Entre las elecciones de norma más populares para el problema de regresión con regularización están:
donde: $I$ es la matriz identidad. Este problema siempre tiene solución (aún si $A$ es de $rank$ incompleto) pues $A^TA + \lambda I$ es una matriz definida positiva para $\lambda >0$. La solución está dada por: $\beta^* = (A^TA + \lambda I)^{-1}A^Ty$.
Comentario: es posible probar que los problemas anteriores son equivalentes a:
$$\displaystyle \min_{\beta \in \mathbb{R}^n} ||A\beta-y||_2^2$$para el caso de ridge y:
para el caso de lasso y con $t$ un parámetro que define la regularización y está relacionado con $\lambda$. Las formulaciones anteriores ayudan a visualizar lo que en el proceso de optimización se está buscando:
en el dibujo anterior las curvas de nivel de la función objetivo (convexa) se representan como elipses y la variable de optimización es $\beta \in \mathbb{R}^2$. Del lado izquierdo tenemos la bola unitaria bajo la norma $1$ que corresponde a la regularización lasso y del lado derecho la bola unitaria bajo la norma $2$ que corresponde a la regularización ridge. En ambos dibujos se observa que la solución está dada por $\beta^*$ y que resulta de la intersección de la curva de nivel que toca a la bola unitaria respectiva.
para valores $\alpha \in [0,1]$. Obsérvese si $\alpha = 0$ se tiene la regularización ridge y si $\alpha=1$ se tiene la regularización lasso. Este tipo de regularización realiza un equilibrio entre ambas regularizaciones.
In [1]:
install.packages(c("latex2exp","glmnet"),lib="/usr/local/lib/R/site-library/",
repos="https://cran.itam.mx/")
En este primer ejemplo no usamos regularización, es un problema de mínimos cuadrados lineales.
In [2]:
#load numerical differentiation
#load utils
#load algorithms for unconstrained convex optimization
#load line search
dir_R="algoritmos/R"
source(paste(dir_R,"/numerical_differentiation.R",
sep=""))
source(paste(dir_R,"/utils.R",
sep=""))
source(paste(dir_R,"/algorithms_for_uco.R",
sep=""))
source(paste(dir_R,"/line_search.R",
sep=""))
In [3]:
library(ggplot2)
library(latex2exp)
library(glmnet)
library(magrittr)
library(dplyr)
Generamos puntos pseudo aleatorios:
In [4]:
set.seed(1989) #para reproducibilidad
mpoints <- 20
df <- data.frame(x=rnorm(mpoints))
y <- -3*df$x + rnorm(mpoints,2,1)
df$y <- y
In [5]:
gg <- ggplot(data=df, aes(x=x, y=y))
In [6]:
gg +
geom_point(aes(x=x,y=y),size=2)
Usamos la función lm del paquete stats de R para ajustar un modelo de regresión lineal:
In [7]:
linear_model <- lm(df$y~df$x)
In [8]:
print(linear_model$coefficients)
In [9]:
gg +
geom_point(aes(x=x,y=y),size=2) +
geom_smooth(method='lm',colour='red')
Aplicamos el método de descenso en gradiente para comparar con lo calculado vía lm
Recordamos que el problema de optimización es: $$\displaystyle \min_{\beta \in \mathbb{R}^n} \quad \frac{1}{2}y^Ty-\beta^TA^Ty + \frac{1}{2}\beta^TA^TA\beta$$
Función objetivo:
In [10]:
cte <- sum(y*y)
In [11]:
A <- matrix(c(rep(1,mpoints),df$x),nrow=mpoints)
In [12]:
fo <-function(beta)1/2*cte - sum(beta*(t(A)%*%y)) + 1/2*sum(beta*(t(A)%*%(A%*%beta)))
#obsérvese que no se realiza el producto A^TA
Punto inicial $\beta^{(0)}:$
In [13]:
beta_0 <- matrix(c(0,0),nrow=2)
In [14]:
beta_ast <- c(linear_model$coefficients[1],linear_model$coefficients[2])
$\beta^*$ (punto óptimo):
In [15]:
print(beta_ast)
$p^*$ (valor óptimo):
In [16]:
p_ast <- fo(beta_ast)
In [17]:
p_ast
argumentos para el método de descenso en gradiente:
In [18]:
tol <- 1e-8
tol_backtracking <- 1e-14
maxiter <- 30
In [19]:
l<-gradient_descent(fo, beta_0, tol, tol_backtracking, beta_ast, p_ast, maxiter)
Soluciones que están en la lista l:
In [20]:
beta <- l[[1]]
total_of_iterations <- l[[2]]
Err_plot <- l[[3]]
beta_plot <- l[[4]]
$\beta$ aproximada por el método de descenso en gradiente:
In [21]:
print(beta)
$\beta^*$:
In [22]:
print(beta_ast)
Error relativo:
In [23]:
compute_error(beta_ast, beta)
Tenemos alrededor de $7$ dígitos de precisión.
Secuencia de minimización $\beta^{(k)}$:
In [24]:
beta_plot
In [25]:
total_of_iterations
Gráfica de error relativo:
In [26]:
gg <- ggplot()
In [27]:
gg +
geom_line(aes(x=25:total_of_iterations,y=Err_plot[25:length(Err_plot)])) +
xlab('Iterations') + ylab(TeX('Error relativo entre f_o(x^k) y p^*'))
En el paquete de R glmnet package se tiene la función del mismo nombre glmnet en la que se ajustan modelos lineales generalizados con penalización elastic net (compromiso entre lasso y ridge). La función objetivo en el caso de regresión lineal (familia Gaussiana) utiliza una pérdida cuadrática con regularización elastic net:
$$\displaystyle \min_{(\beta_0, \beta) \in \mathbb{R}^{n+1}} \frac{1}{2m} \sum_{i=1}^m(y_i -\beta_0 - x_i^T\beta)^2 + \lambda \left (\frac{(1-\alpha)}{2}||\beta||^2_2 + \alpha ||\beta||_1 \right )$$donde: $x_i \in \mathbb{R}^n$.
Véase el artículo Regularization Paths for Generalized Linear Models via Coordinate Descent para esta formulación.
Obsérvese que no se penaliza la variable $\beta_0$.
Comentarios:
Lo que continúa en la nota son comparaciones entre lo obtenido por el paquete de glmnet de R vs implementaciones simples de los métodos de descenso por coordenadas y Newton realizadas por el prof. La implementación realizada en el paquete es mucho más general y eficiente que lo realizado por el prof. No se pretende realizar comparaciones en tiempo, memoria ni generalidad en la solución de problemas. Aún así, lo presentado a continuación ayuda a entender el problema que se resuelve y la metodología utilizada.
También en los ejemplos no se realizará estandarización de variables (aunque es recomendable realizar esto para ejemplos reales...).
En este segundo ejemplo utilizamos la regularización lasso.
Obsérvese que para este caso la función objetivo en glmnet es de la forma:
Comentario: recuérdese que $||\beta||_1 = \displaystyle \sum_{i=1}^n |\beta_i|$ por lo que la función objetivo continúa siendo convexa pero no es diferenciable en el vector $\beta = 0$ (el valor absoluto es una función no diferenciable en el punto $0$).
Simulamos algunos datos:
Obs: se utilizarán los siguientes argumentos para la función de R:
In [28]:
reg<-.5 #este es lambda * alpha, el parámetro de regularización
#en la formulación de glmnet
In [29]:
fit <- glmnet(A,y,alpha=1,lambda=reg,standardize=F,nlambda=1,thresh=1e-8)
In [30]:
beta_ast <- as.matrix(fit$beta)
$\beta^*$:
In [31]:
print(beta_ast)
Obsérvese que $\beta^*_0$ es $0$ y por tanto se puede eliminar el intercepto del modelo.
Comentario: A continuación se elimina la primera columna de la matriz $A$ por la observación anterior. La primer columna de $A$ (columna de $1$s) implica considerar un modelo con intercepto: $\beta_0$. Además como se mencionó anteriormente, la implementación en glmnet para el caso de lasso es más general y tiende a hacer $0$s los coeficientes estimados. La implementación del prof no está realizando esto (pues el objetivo es mostrar el uso de métodos de descenso para resolver diferentes problemas, el objetivo no es obtener los mismos resultados que glmnet) por lo que los coeficientes que son $0$ no serán estimados correctamente en esta implementación.
In [32]:
print(head(A))
Eliminamos la primer columna de $A$ que ayuda a la estimación del intercepto:
In [33]:
A<-A[,-1]
In [34]:
print(head(A))
El objetivo entonces es estimar una sola $\beta$: $\beta_1$.
In [35]:
beta_ast<-beta_ast[2]
In [36]:
print(beta_ast)
Usaremos el método de descenso por coordenadas y el método de Newton para aproximar a $\beta_1^*$. Ver 4.2.Descenso_por_coordenadas_R.ipynb, 4.2.Metodo_de_Newton_Python. Además usaremos la siguiente función quita_signo que ayuda a aproximar la derivada de la función objetivo $f_o$ y así lidiar con la no diferenciabilidad en $0$.
La función quita_signo realiza:
obsérvese que el valor elegido en el tercer caso es el más pequeño positivo normalizado en un sistema de punto flotante, ver 1.2.Sistema_de_punto_flotante.
In [37]:
quita_signo<-function(beta){
beta<-sign(beta)*beta
#la siguiente variable es un índice que localiza aquellas entradas del
#vector beta en valor absoluto que son cercanas a 0.
ind <- beta < .Machine$double.xmin & beta > -.Machine$double.xmin
#se asigna a cada entrada localizada en ind el valor más pequeño normalizado
#en un sistema de punto flotante
beta[ind] <- .Machine$double.xmin
beta
}
In [38]:
.Machine$double.xmin
Comentario: la definición de la función quita_signo se basa en lo que se conoce como subdifferential que es un conjunto de subderivatives, útiles para generalizar las derivadas para funciones convexas que no son diferenciables en puntos de su dominio.
Así, la función objetivo es:
In [39]:
fo <-function(beta)1/mpoints*(1/2*cte - sum(beta*(A*y)) +
1/2*sum(beta*(A*(A*beta)))) +
reg*sum(quita_signo(beta))
Valor óptimo:
In [40]:
p_ast <- fo(beta_ast)
In [41]:
p_ast
Punto inicial $\beta^{(0)}:$
In [42]:
beta_0<-0
Argumentos para el método de descenso por coordenadas:
In [43]:
tol <- 1e-8
tol_backtracking <- 1e-14
maxiter <- 30
In [44]:
l<-coordinate_descent(fo, beta_0, tol, tol_backtracking, beta_ast, p_ast, maxiter)
Soluciones que están en la lista l:
In [45]:
beta <- l[[1]]
total_of_iterations <- l[[2]]
Err_plot <- l[[3]]
beta_plot <- l[[4]]
$\beta$ aproximada por el método de descenso por coordenadas:
In [46]:
print(beta)
$\beta^*$:
In [47]:
print(beta_ast)
Error relativo:
In [48]:
compute_error(beta_ast, beta)
Tenemos alrededor de $2$ dígitos de precisión.
Secuencia de minimización $\beta^{(k)}$:
In [49]:
print(beta_plot)
In [50]:
gg +
geom_point(aes(x=beta_plot,y=0),size=2) +
annotate(geom='text', x=0, y=0,
label=TeX("x^{(0)}", output='character'), parse=TRUE) +
xlab('x') + ylab('y') +
ggtitle(TeX('Iter del método de descenso por coordenadas para $f_o$'))
In [51]:
l<-Newtons_method(fo, beta_0, tol, tol_backtracking, beta_ast, p_ast, maxiter)
Soluciones que están en la lista l:
In [52]:
beta <- l[[1]]
total_of_iterations <- l[[2]]
Err_plot <- l[[3]]
beta_plot <- l[[4]]
$\beta$ aproximada por el método de Newton:
In [53]:
print(beta)
$\beta^*$:
In [54]:
print(beta_ast)
Error relativo:
In [55]:
compute_error(beta_ast, beta)
Tenemos alrededor de $2$ dígitos de precisión.
Secuencia de minimización $\beta^{(k)}$:
In [56]:
print(beta_plot)
In [57]:
gg +
geom_point(aes(x=beta_plot,y=0),size=2) +
annotate(geom='text', x=0, y=0,
label=TeX("x^{(0)}", output='character'), parse=TRUE) +
xlab('x') + ylab('y') +
ggtitle(TeX('Iter del método de Newton para $f_o$'))
Comentario: En ambos métodos se aproxima de forma correcta a $\beta_1^*$.
Simulamos otros datos:
In [58]:
set.seed(1989) #para reproducibilidad
mpoints <- 50
x1 <- rnorm(mpoints)
x2 <- rnorm(mpoints,2,1)
y <- 3*x1 -.5*x2
In [59]:
A<-cbind(x1,x2)
In [60]:
print(head(A))
In [61]:
print(head(y))
Reconstruímos a la función objetivo con el nuevo valor de la constante y el número de puntos:
In [62]:
cte <- sum(y*y)
mpoints<-nrow(A)
In [63]:
cte
In [64]:
mpoints
In [65]:
fo <-function(beta)1/mpoints*(1/2*cte - sum(beta*(t(A)%*%y)) +
1/2*sum(beta*(t(A)%*%(A%*%beta)))) +
reg*sum(quita_signo(beta))
Solución vía glmnet sin intercepto:
In [66]:
fit <- glmnet(A,y,alpha=1,lambda=reg,standardize=F,nlambda=1,intercept=F,thresh=1e-8)
In [67]:
beta_ast <- as.matrix(fit$beta)
In [68]:
print(beta_ast)
Solución vía método de Newton:
In [69]:
beta_0<-c(1,1)
In [70]:
tol <- 1e-8
tol_backtracking <- 1e-14
maxiter <- 30
p_ast <- fo(beta_ast)
In [71]:
l<-Newtons_method(fo, beta_0, tol, tol_backtracking, beta_ast, p_ast, maxiter)
In [72]:
beta <- l[[1]]
total_of_iterations <- l[[2]]
Err_plot <- l[[3]]
beta_plot <- l[[4]]
$\beta$ aproximada por el método de Newton:
In [73]:
print(beta)
$\beta^*$:
In [74]:
print(beta_ast)
Error relativo:
In [75]:
compute_error(beta_ast, beta)
In [76]:
print(head(mtcars))
Utilizamos las variables numéricas disp y drat:
In [77]:
y <- mtcars %>% select(mpg) %>% as.matrix()
X <- mtcars %>% select(-mpg) %>% as.matrix()
In [78]:
A<-X[,c(2,4)]
In [79]:
print(head(A))
Ajuste vía glmnet:
In [80]:
fit <- glmnet(A,y,alpha=1,lambda=reg,standardize=F,nlambda=1,intercept=F,thresh=1e-8)
In [81]:
beta_ast <- as.matrix(fit$beta)
$\beta^*$:
In [82]:
print(beta_ast)
Solución vía método de Newton
Reconstruímos a la función objetivo con el nuevo valor de la constante y el número de puntos:
In [83]:
cte <- sum(y*y)
mpoints<-nrow(A)
In [84]:
cte
In [85]:
mpoints
In [86]:
beta_0<-c(1,1)
In [87]:
tol <- 1e-8
tol_backtracking <- 1e-14
maxiter <- 30
p_ast <- fo(beta_ast)
In [88]:
l<-Newtons_method(fo, beta_0, tol, tol_backtracking, beta_ast, p_ast, maxiter)
In [89]:
beta <- l[[1]]
total_of_iterations <- l[[2]]
Err_plot <- l[[3]]
beta_plot <- l[[4]]
$\beta$ aproximada por el método de Newton:
In [90]:
print(beta)
$\beta^*$:
In [91]:
print(beta_ast)
Error relativo:
In [92]:
compute_error(beta_ast, beta)
Tenemos alrededor de $5$ dígitos de precisión.
In [93]:
reg<-0.2
Ajuste vía glmnet:
In [94]:
fit <- glmnet(A,y,alpha=1,lambda=reg,standardize=F,nlambda=1,intercept=F,thresh=1e-8)
$\beta^*:$
In [95]:
beta_ast <- as.matrix(fit$beta)
In [96]:
print(beta_ast)
Solución vía método de Newton:
In [97]:
beta_0<-c(1,1)
In [98]:
tol <- 1e-8
tol_backtracking <- 1e-14
maxiter <- 30
p_ast <- fo(beta_ast)
In [99]:
l<-Newtons_method(fo, beta_0, tol, tol_backtracking, beta_ast, p_ast, maxiter)
In [100]:
beta <- l[[1]]
total_of_iterations <- l[[2]]
Err_plot <- l[[3]]
beta_plot <- l[[4]]
$\beta$ aproximada por el método de Newton:
In [101]:
print(beta)
$\beta^*:$
In [102]:
print(beta_ast)
Error relativo:
In [103]:
compute_error(beta_ast, beta)
Tenemos alrededor de $5$ dígitos de precisión.
Solución vía descenso por coordenadas:
In [104]:
beta_0<-c(0,0)
In [105]:
l<-coordinate_descent(fo, beta_0, tol, tol_backtracking, beta_ast, p_ast, maxiter)
In [106]:
beta <- l[[1]]
total_of_iterations <- l[[2]]
Err_plot <- l[[3]]
beta_plot <- l[[4]]
$\beta$ aproximada por el método de descenso por coordenadas:
In [107]:
print(beta)
$\beta^*:$
In [108]:
print(beta_ast)
Error relativo:
In [109]:
compute_error(beta_ast, beta)
Tenemos alrededor de $2$ dígitos de precisión.
Secuencia de minimización:
In [110]:
beta_plot
In [111]:
gg +
geom_point(aes(x=beta_plot[1,],y=beta_plot[2,]),size=2) +
annotate(geom='text', x=0, y=-0.1,
label=TeX("x^{(0)}", output='character'), parse=TRUE) +
xlab('x') + ylab('y') +
ggtitle(TeX('Iter del método de descenso por coordenadas para $f_o$'))
In [112]:
l<-gradient_descent(fo, beta_0, tol, tol_backtracking, beta_ast, p_ast, maxiter)
In [113]:
beta <- l[[1]]
total_of_iterations <- l[[2]]
Err_plot <- l[[3]]
beta_plot <- l[[4]]
$\beta$ aproximada por el método de descenso en gradiente:
In [114]:
print(beta)
$\beta^*:$
In [115]:
print(beta_ast)
Error relativo:
In [116]:
compute_error(beta_ast,beta)
Tenemos un error del $97\%$!.
In [117]:
beta_plot
In [118]:
gg +
geom_point(aes(x=beta_plot[1,],y=beta_plot[2,]),size=2) +
annotate(geom='text', x=0, y=-0.01,
label=TeX("x^{(0)}", output='character'), parse=TRUE) +
xlab('x') + ylab('y') +
ggtitle(TeX('Iter del método de descenso en gradiente para $f_o$'))
In [119]:
A<-X[,c(2,4,5,6)]
In [120]:
print(head(A))
Solución vía glmnet:
In [121]:
fit <- glmnet(A,y,alpha=1,lambda=reg,standardize=F,nlambda=1,intercept=F,thresh=1e-8)
$\beta^*$:
In [122]:
beta_ast <- as.matrix(fit$beta)
In [123]:
print(beta_ast)
Solución vía método de Newton:
In [124]:
beta_0<-c(1,1,1,1)
In [125]:
tol <- 1e-8
tol_backtracking <- 1e-14
maxiter <- 30
p_ast <- fo(beta_ast)
Newtons method
In [126]:
l<-Newtons_method(fo, beta_0, tol, tol_backtracking, beta_ast, p_ast, maxiter)
In [127]:
beta <- l[[1]]
total_of_iterations <- l[[2]]
Err_plot <- l[[3]]
beta_plot <- l[[4]]
$\beta$ aproximada por el método de Newton:
In [128]:
print(beta)
$\beta^*$:
In [129]:
print(beta_ast)
Error relativo:
In [130]:
compute_error(beta_ast, beta)
Tenemos alrededor de $3$ dígitos de precisión.
Secuencia de minimización:
In [131]:
beta_plot
Función objetivo en glmnet:
Obsérvese que para este caso la función objetivo en glmnet es de la forma ($\alpha=0$):
Comentarios:
A diferencia del caso de lasso, la función objetivo para la regularización ridge sí es diferenciable en cualquier punto de su dominio. Recuérdese que $||\beta||^2_2 = \displaystyle \sum_{i=1}^n \beta_i^2$ y por tanto la función objetivo continúa siendo convexa.
También en este caso a diferencia de lasso la solución está dada por la expresión vía la descomposición en valores singulares de la matriz $A$, ver 3.3.d.SVD.ipynb:
donde: $D$ es una matriz diagonal con entradas $\frac{\sigma_i}{\sigma_i^2 + \lambda}$ para $i=1,\dots , n$.
Usamos el siguiente parámetro de regularización:
In [132]:
reg<-.5 #este es lambda, el parámetro de regularización
#en la formulación de glmnet
In [133]:
y <- mtcars %>% select(mpg) %>% as.matrix()
X <- mtcars %>% select(-mpg) %>% as.matrix()
In [134]:
A<-X[,c(2,4)]
In [135]:
print(head(A))
Solución vía glmnet:
In [136]:
fit <- glmnet(A,y,alpha=0,lambda=reg,standardize=F,nlambda=1,intercept=F,thresh=1e-8)
In [137]:
beta_ast <- as.matrix(fit$beta)
$\beta^*:$
In [138]:
print(beta_ast)
Solución vía SVD:
In [139]:
#svd of A
singular_value_decomposition <- svd(A)
s <- singular_value_decomposition$d
u <- singular_value_decomposition$u
v <- singular_value_decomposition$v
cte_svd <- s/(s^2+reg)*(t(u)%*%y)
In [140]:
beta_ridge <- v%*%cte_svd
In [141]:
print(beta_ridge)
In [142]:
compute_error(beta_ast,beta_ridge)
Tenemos aproximadamente dos dígitos de precisión.
In [143]:
A<-X[,c(2,4,5,6)]
In [144]:
print(head(A))
Solución vía glmnet:
In [145]:
fit <- glmnet(A,y,alpha=0,lambda=reg,standardize=F,nlambda=1,intercept=F,thresh=1e-8)
In [146]:
beta_ast <- as.matrix(fit$beta)
$\beta^*$:
In [147]:
print(beta_ast)
Solución vía SVD:
In [148]:
#svd of A
singular_value_decomposition <- svd(A)
s <- singular_value_decomposition$d
u <- singular_value_decomposition$u
v <- singular_value_decomposition$v
cte_svd <- s/(s^2+reg)*(t(u)%*%y)
In [149]:
beta_ridge <- v%*%cte_svd
In [150]:
print(beta_ridge)
Error relativo:
In [151]:
compute_error(beta_ast,beta_ridge)
Tenemos $3\%$ de error. La razón de la mala estimación es debido al mal condicionamiento de la matriz $A^TA$:
In [152]:
kappa(t(A)%*%A,exact=TRUE)
El mal condicionamiento de $A$ en Estadística lo relacionan con problemas de multicolinealidad, véase Multicollinearity:
In [153]:
cor(A)
Entonces tenemos alta correlación entre la variable disp y wt.
In [154]:
A<-X[,c(4,5,6)]
In [155]:
kappa(t(A)%*%A, exact=TRUE)
In [156]:
cor(A)
Solución vía glmnet:
In [157]:
fit <- glmnet(A,y,alpha=0,lambda=reg,standardize=F,nlambda=1,intercept=F,thresh=1e-8)
In [158]:
beta_ast <- as.matrix(fit$beta)
$\beta^*$:
In [159]:
print(beta_ast)
Solución vía SVD:
In [160]:
#svd of A
singular_value_decomposition <- svd(A)
s <- singular_value_decomposition$d
u <- singular_value_decomposition$u
v <- singular_value_decomposition$v
cte_svd <- s/(s^2+reg)*(t(u)%*%y)
In [161]:
beta_ridge <- v%*%cte_svd
In [162]:
print(beta_ridge)
Error relativo:
In [163]:
compute_error(beta_ast,beta_ridge)
Tenemos alrededor de $2$ dígitos de precisión.
Referencias:
Regularization Paths for Generalized Linear Models via Coordinate Descent
Ver 3_minimos_cuadrados para una introducción al problema de mínimos cuadrados con ejemplos en Python3.