====== 10 - Regularización ====== La regularización consiste en añadir una penalización a la función de coste para minimizar los pesos con el objeto de hacer el modelo más sencillo. Las regularizaciones más usadas en ML son: * Ridge (L2)\\ \\ * Lasso (L1)\\ \\ * ElasticNet\\ \\ Para ver como funcionan, vamos a crear un modelo de regresión lineal con los datos del dataset {{ :clase:ia:saa:6_optimizacion_modelo:meatspec.zip |}}. Este dataset contiene datos del contenido de grasa en carnes según 100 parámetros. El objetivo es predecir ese contenido de grasa. Librerías necesarias: import pandas as pd import numpy as np from sklearn.model_selection import train_test_split from sklearn.linear_model import LinearRegression from sklearn.metrics import mean_squared_error from sklearn.preprocessing import StandardScaler import matplotlib.pyplot as plt from sklearn.linear_model import Ridge from sklearn.linear_model import Lasso from sklearn.linear_model import ElasticNet #Regularización con Cross-Validate para encontrar el mejor valor de alpha from sklearn.linear_model import RidgeCV from sklearn.linear_model import LassoCV from sklearn.linear_model import ElasticNetCV # Configuración warnings # ============================================================================== import warnings warnings.filterwarnings('ignore') Cargamos los datos del dataset: df = pd.read_csv('./../../../Datasets/meatspec.csv', index_col=0) #Renombramos la columna objetivo (fat) a target df.rename(columns = {'fat': 'target'}, inplace = True) Separamos entrenamiento y test: X = df.drop(columns = ['target'], axis = 1) y = df['target'] standard_scaler = StandardScaler() X_scaled = standard_scaler.fit_transform(X) X_scaled = pd.DataFrame(X_scaled, columns = X.columns) X_train, X_test, y_train, y_test = train_test_split(X_scaled, y, random_state=9) Para que funcione bien la regularización los datos tienen que estar normalizados. Aunque deberíamos normalizar sólo los datos de train y crear una pipeline para normalizar los de test, lo hacemos así para simplificar. Creamos el modelo: model_lr = LinearRegression() _ = model_lr.fit(X_train, y_train) Miramos el error (en este caso RMSE): #Error de entrenamiento y_predict = model_lr.predict(X_train) train_error = mean_squared_error(y_train, y_predict, squared = False) print("RMSE entrenamiento: %.2f" % train_error) #Error de test y_predict = model_lr.predict(X_test) ols_error = mean_squared_error(y_test, y_predict, squared = False) print("RMSE test: %.2f" % ols_error) RMSE entrenamiento: 0.77 RMSE test: 3.69 Como vemos, el error de entrenamiento se aleja bastante del error de test, síntoma de sobreajuste del modelo. Vamos a mostrar un gráfico para ver los coeficientes del modelo. Para eso, creamos una función, ya que la utilizaremos más adelante: def plot_coefficients(model, X): df_coefficients = pd.DataFrame( { 'predictor': X.columns, 'coef': model.coef_ } ) figure = plt.figure(figsize=(15, 5)) axes = figure.add_subplot() axes.stem(df_coefficients['predictor'], df_coefficients['coef'], markerfmt=' ') axes.tick_params(axis = 'x', labelrotation = 90, labelsize = 7) axes.set_xlabel('variable') axes.set_ylabel('coeficientes') axes.set_title('Coeficientes del modelo') plot_coefficients(model_lr, X_train) {{ :clase:ia:saa:6_optimizacion_modelo:coef_1.png?800 |}} El rango de nuestros coeficientes es [-30000, 30000]. ===== Ridge (L2) ===== La regularización Ridge penaliza la suma de los coeficientes elevados al cuadrado: $$ Error=\frac 1 n \displaystyle\sum_{i=1}^{n} {(y_i-y'_i)^2 } + \alpha \frac 1 {m} \displaystyle\sum_{j=1}^{m} {|w_j|^2} $$ A esta penalización se le conoce como l2 y tiene el efecto de reducir de forma proporcional el valor de todos los coeficientes del modelo pero sin que estos lleguen a cero. El grado de penalización está controlado por el hiperparámetro α. Cuando α = 0 , la penalización es nula y el resultado es equivalente al de un modelo lineal por mínimos cuadrados ordinarios (OLS). A medida que α aumenta, mayor es la penalización y menor el valor de los predictores. Para aplicar Ridge, vamos a utilizar la función de //sklearn// [[https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Ridge.html|Ridge]]. Primero, vamos a comprobar que con α = 0 el modelo se queda igual: model_ridge = Ridge(alpha=0) _ = model_ridge.fit(X_train, y_train) plot_coefficients(model_ridge, X_train) {{ :clase:ia:saa:6_optimizacion_modelo:coef_1.png?600 |}} Como vemos en el gráfico, los coeficientes no han variado. Vamos a probar con α = 0.5: model_ridge = Ridge(alpha=0.5) _ = model_ridge.fit(X_train, y_train) plot_coefficients(model_ridge, X_train) {{ :clase:ia:saa:6_optimizacion_modelo:coef_2.png?600 |}} Ahora sí que ha bajado mucho el rango de los coeficientes, [-6, 10] más o menos. Vemos el error: #Error de entrenamiento y_predict = model_ridge.predict(X_train) train_error = mean_squared_error(y_train, y_predict, squared = False) print("RMSE entrenamiento: %.2f" % train_error) #Error de test y_predict = model_ridge.predict(X_test) test_error = mean_squared_error(y_test, y_predict, squared = False) print("RMSE test: %.2f" % test_error) RMSE entrenamiento: 3.27 RMSE test: 2.68 Nuestro error de entrenamiento ha empeorado, pero el de test ha mejorado. Además, ya no hay tanta diferencia entre ambos, con lo que hemos reducido el sobreajuste. ¿Y con α = 1? model_ridge = Ridge(alpha=1) _ = model_ridge.fit(X_train, y_train) plot_coefficients(model_ridge, X_train) {{ :clase:ia:saa:6_optimizacion_modelo:coef_3.png?600 |}} #Error de entrenamiento y_predict = model_ridge.predict(X_train) train_error = mean_squared_error(y_train, y_predict, squared = False) print("RMSE entrenamiento: %.2f" % train_error) #Error de test y_predict = model_ridge.predict(X_test) test_error = mean_squared_error(y_test, y_predict, squared = False) print("RMSE test: %.2f" % test_error) RMSE entrenamiento: 3.52 RMSE test: 2.88 El valor de los coeficientes baja todavía más, pero el error de test también ha aumentado. ¿Qué valor de α es el adecuado? Al ser un hiperparámetro que tenemos que definir a la hora crear el modelo, necesitamos saber cual es el valor óptimo de α. Para eso, podemos utilizar la función [[https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.RidgeCV.html|RidgeCV]], que buscará mediante validación cruzada el valor óptimo de α que minimice la función de coste (en nuestro caso RMSE). Esta función almacena en el atributo *alpha_* el mejor valor de α encontrado. Además, si ponemos el parámetro *store_cv_values = True* almacenará los errores para cada valor de α probado. #np.logspace devuelve números espaciados uniformemente en una escala logarítmica. model_ridgeCV = RidgeCV( alphas = np.logspace(-10, 2, 200), store_cv_values = True ) _ = model_ridgeCV.fit(X_train, y_train) Creamos una función para ver los coeficientes del modelo en función de la regularización: def plot_coefficients_reg(alphas, model, X, y): coefs = [] for alpha in alphas: if model == 'ridge': model_temp = Ridge(alpha=alpha) elif model == 'lasso': model_temp = Lasso(alpha=alpha) elif model == 'elastic': model_temp = ElasticNet(alpha=alpha) model_temp.fit(X, y) coefs.append(model_temp.coef_) figure = plt.figure(figsize=(15, 5)) axes = figure.add_subplot() axes.plot(alphas, coefs) axes.set_xscale('log') axes.set_xlabel('alpha') axes.set_ylabel('coeficientes') axes.set_title('Coeficientes del modelo en función de la regularización') plot_coefficients_reg(model_ridgeCV.alphas, 'ridge', X_train, y_train) {{ :clase:ia:saa:6_optimizacion_modelo:coef_reg_1.png?600 |}} Como vemos, cuanto más alto es el valor de α, más se acercan los coeficientes a 0. Por último, creamos otra función para ver la evolución del error en función del valor de α: def plot_error(cv_values, alphas): # np.sqrt para trabajar con RMSE # axis = 0 para calcular la media de cada columna mse_mean = np.sqrt(cv_values.mean(axis = 0)) #Sacamos el óptimo best_mean = np.min(mse_mean) #index_best devuelve un array de arrays index_best = np.where(mse_mean == best_mean)[0][0] best_alpha = alphas[index_best] figure = plt.figure(figsize=(15, 5)) axes = figure.add_subplot() axes.plot(alphas, mse_mean) axes.set_xscale('log') axes.set_ylim([0,None]) axes.set_title('Evolución del error CV en función de la regularización') axes.set_xlabel('alpha') axes.set_ylabel('RMSE') label = 'best alpha: ' + str(best_alpha) axes.axvline(x = best_alpha, c = "blue", linestyle = '--', label = label) axes.legend(fontsize=15,facecolor='#CDCDCD',labelcolor="#000000") # modelo.cv_values almacena el mse de cv para cada valor de alpha. Tiene # dimensiones (n_samples, n_alphas) plot_error(model_ridgeCV.cv_values_, model_ridgeCV.alphas) {{ :clase:ia:saa:6_optimizacion_modelo:error1.png?600 |}} Mostramos el mejor valor de α encontrado por la función y el valor de los coeficientes: print(f"El mejor valor de alpha encontrado es {model_ridgeCV.alpha_}") El mejor valor de alpha encontrado es 2.3272024789604072e-05 plot_coefficients(model_ridgeCV, X_train) {{ :clase:ia:saa:6_optimizacion_modelo:coef_4.png?600 |}} Error del modelo: #Error de entrenamiento y_predict = model_ridgeCV.predict(X_train) train_error = mean_squared_error(y_train, y_predict, squared = False) print("RMSE entrenamiento: %.2f" % train_error) #Error de test y_predict = model_ridgeCV.predict(X_test) ridge_error = mean_squared_error(y_test, y_predict, squared = False) print("RMSE test: %.2f" % ridge_error) RMSE entrenamiento: 1.91 RMSE test: 1.80 Nuestros errores se han acercado y además son más bajos que en la regresión lineal. Además, los coeficientes finales tienen un rango mucho más bajo ([-850, 650] más o menos). Aunque el valor óptimo de α es aquel con el que se minimiza el error de validación cruzada, una práctica extendida es utilizar, en lugar de este, el mayor valor de α que se aleja menos de una desviación típica del óptimo. De este modo, se consigue un modelo más sencillo, pero cuya capacidad predictiva es similar a la conseguida con el modelo más complejo. ===== Lasso (L1) ===== Al contrario que Ridge, Lasso minimiza los coeficientes pudiendo llegar a 0, con lo que puede que alguno de ellos quede eliminado del modelo. $$ Error=\frac 1 n \displaystyle\sum_{i=1}^{n} {(y_i-y'_i)^2 } + \alpha \frac 1 m \displaystyle\sum_{j=1}^{m} {|w_j|^1} $$ Igual que hemos hecho con Ridge, vamos a utilizar [[http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LassoCV.html|LassoCV]] para encontrar el mejor valor de α: model_lassoCV = LassoCV( alphas = np.logspace(-10, 3, 200), cv = 10 ) _ = model_lassoCV.fit(X_train, y_train) Vemos que si aumentamos el valor de α, los coeficientes se acercan a 0 (al contrario que con Ridge, aquí sí que pueden llegar a 0, con lo que es una forma de reducir la dimensionalidad de nuestro modelo). plot_coefficients_reg(model_lassoCV.alphas_, 'lasso', X_train, y_train) {{ :clase:ia:saa:6_optimizacion_modelo:plot_coeficients1.png?600 |}} Mostramos la variación del error en función de α y el mejor valor encontrado: plot_error(model_lassoCV.mse_path_.T, model_lassoCV.alphas_) {{ :clase:ia:saa:6_optimizacion_modelo:plot_errors1.png?600 |}} print(f"Mejor valor de alpha encontrado: {model_lassoCV.alpha_}") Mejor valor de alpha encontrado: 1e-10 Miramos el error de entrenamiento y test: #Error de entrenamiento y_predict = model_lassoCV.predict(X_train) train_error = mean_squared_error(y_train, y_predict, squared = False) print("RMSE entrenamiento: %.2f" % train_error) #Error de test y_predict = model_lassoCV.predict(X_test) lasso_error = mean_squared_error(y_test, y_predict, squared = False) print("RMSE test: %.2f" % lasso_error) RMSE entrenamiento: 4.87 RMSE test: 3.86 Aunque hemos reducido el sobreajuste, no mejoramos casi nada el modelo sin regularizar. Vamos a ver como quedan nuestros coeficientes: # Coeficientes del modelo # ============================================================================== df_coeficientes = pd.DataFrame( {'predictor': X_train.columns, 'coef': model_lassoCV.coef_.flatten()} ) # Predictores incluidos en el modelo (coeficiente != 0) df_coeficientes[df_coeficientes.coef != 0] predictor coef 0 V1 59.484875 1 V2 -16.613824 2 V3 -15.894723 3 V4 -15.250222 4 V5 -14.543708 ... ... ... 95 V96 1.195200 96 V97 -0.134883 97 V98 -1.917606 98 V99 -3.943738 99 V100 -5.612002 100 rows × 2 columns No ha eliminado ningún coeficiente (aunque el valor es más bajo). Vamos a ver el número de coeficientes incluidos en el modelo con Lasso en función del valor de α: # Número de predictores incluidos (coeficiente !=0) en función de alpha # ============================================================================== alphas = model_lassoCV.alphas_ n_predictores = [] for alpha in alphas: modelo_temp = Lasso(alpha=alpha) modelo_temp.fit(X_train, y_train) coef_no_cero = np.sum(modelo_temp.coef_ != 0) n_predictores.append(coef_no_cero) fig, ax = plt.subplots(figsize=(7, 3.84)) ax.plot(alphas, n_predictores) ax.set_xscale('log') ax.set_ylim([-15,None]) ax.set_xlabel('alpha') ax.set_ylabel('nº predictores') ax.set_title('Predictores incluidos en función de la regularización') {{ :clase:ia:saa:6_optimizacion_modelo:coef_alpha.png?600 |}} Como se puede apreciar en la anterior imagen, Lasso empieza a eliminar coeficientes a partir de $\alpha = 10^{-3}$ más o menos. Antes hemos dicho que es habitual elegir el mayor valor de α que se aleja menos de una desviación típica del óptimo. Vamos a comprobar si mejoramos nuestro modelo: # Evolución del error en función de alpha # ============================================================================== # modelo.mse_path_ almacena el mse de cv para cada valor de alpha. Tiene # dimensiones (n_alphas, n_folds) mse_cv = model_lassoCV.mse_path_.mean(axis=1) mse_sd = model_lassoCV.mse_path_.std(axis=1) # Se aplica la raíz cuadrada para pasar de mse a rmse rmse_cv = np.sqrt(mse_cv) rmse_sd = np.sqrt(mse_sd) # Se identifica el óptimo y el óptimo + 1std min_rmse = np.min(rmse_cv) sd_min_rmse = rmse_sd[np.argmin(rmse_cv)] min_rsme_1sd = np.max(rmse_cv[rmse_cv <= min_rmse + sd_min_rmse]) optimo = model_lassoCV.alphas_[np.argmin(rmse_cv)] # Mejor valor alpha encontrado + 1sd # ============================================================================== min_rmse = np.min(rmse_cv) sd_min_rmse = rmse_sd[np.argmin(rmse_cv)] min_rsme_1sd = np.max(rmse_cv[rmse_cv <= min_rmse + sd_min_rmse]) optimo = model_lassoCV.alphas_[np.argmin(rmse_cv)] optimo_1sd = model_lassoCV.alphas_[rmse_cv == min_rsme_1sd] print(f"Mejor valor de alpha encontrado + 1 desviación estándar: {optimo_1sd}") Mejor valor de alpha encontrado + 1 desviación estándar: [0.18896523] Entrenamos nuestro modelo con ese valor de $\alpha$: # Mejor modelo alpha óptimo + 1sd # ============================================================================== model_lasso_1std = Lasso(alpha=optimo_1sd) model_lasso_1std.fit(X_train, y_train) # Coeficientes del modelo # ============================================================================== df_coeficientes = pd.DataFrame( {'predictor': X_train.columns, 'coef': model_lasso_1std.coef_.flatten()} ) # Predictores incluidos en el modelo (coeficiente != 0) df_coeficientes[df_coeficientes.coef != 0] predictor coef 4 V5 -1.002796 5 V6 -8.277915 6 V7 -7.335488 7 V8 -6.294674 8 V9 -5.383142 9 V10 -4.311487 10 V11 -3.376037 11 V12 -2.394135 12 V13 -1.300096 13 V14 -0.081426 37 V38 18.159401 38 V39 20.662721 39 V40 13.594109 40 V41 2.898530 51 V52 -2.733636 52 V53 -5.269952 53 V54 -2.028877 Ahora sí que ha eliminado características, quedándose sólo con 17 (ha eliminado 83). plot_coefficients(model_lasso_1std, X_train) {{ :clase:ia:saa:6_optimizacion_modelo:coef_5.png?600 |}} Error del modelo: #Error de entrenamiento y_predict = model_lasso_1std.predict(X_train) train_error= mean_squared_error(y_train, y_predict, squared = False) print("RMSE entrenamiento: %.2f" % train_error) #Error de test y_predict = model_lasso_1std.predict(X_test) lasso_1std_error = mean_squared_error(y_test, y_predict, squared = False) print("RMSE test: %.2f" % lasso_error) RMSE entrenamiento: 4.77 RMSE test: 3.86 Aunque nuestro error ha empeorado, hay que tener en cuenta que el modelo se ha simplificado bastante, al eliminar 83 características. ===== ElasticNet ===== Es la unión de Lasso y Ridge. El hiperparámetro **l1_ratio** indica que importancia le da a cada regularización: * 0: Sólo se hace regularización Ridge * 1: Sólo se hace regularización Lasso * 0,5: Se le da la misma importancia a ambas regularizaciones Usamos (como antes) [[https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.ElasticNetCV.html|ElasticNetCV]], pero esta vez tendremos que encontrar dos parámetros óptimos: $\alpha$ y //l1_ratio//: # Creación y entrenamiento del modelo (con búsqueda por CV del valor óptimo alpha) # ============================================================================== # Por defecto ElasticNetCV utiliza el mean squared error model_elasticCV = ElasticNetCV( l1_ratio = [0, 0.1, 0.5, 0.7, 0.9, 0.95, 0.99], alphas = np.logspace(-10, 3, 200), cv = 10 ) _ = model_elasticCV.fit(X = X_train, y = y_train) print(f"Mejor alfa encontrado: {model_elasticCV.alpha_}") print(f"Mejor l1 ratio encontrado: {model_elasticCV.l1_ratio_}") Mejor alfa encontrado: 1e-10 Mejor l1 ratio encontrado: 0.7 Coeficientes del modelo: plot_coefficients(model_elasticCV, X_train) {{ :clase:ia:saa:6_optimizacion_modelo:coef_6.png?600 |}} # Coeficientes del modelo # ============================================================================== df_coeficientes = pd.DataFrame( {'predictor': X_train.columns, 'coef': model_elasticCV.coef_.flatten()} ) # Predictores incluidos en el modelo (coeficiente != 0) df_coeficientes[df_coeficientes.coef != 0] predictor coef 0 V1 59.484876 1 V2 -16.613824 2 V3 -15.894723 3 V4 -15.250222 4 V5 -14.543708 ... ... ... 95 V96 1.195200 96 V97 -0.134883 97 V98 -1.917606 98 V99 -3.943738 99 V100 -5.612002 100 rows × 2 columns En este caso no ha eliminado ningún coeficiente. Podríamos hacer como antes y buscar no los valores óptimos, si no aquellos que simplifiquen el modelo sin que el error se vea muy afectado. Error del modelo: #Error de entrenamiento y_predict = model_elasticCV.predict(X_train) train_error = mean_squared_error(y_train, y_predict, squared = False) print("RMSE entrenamiento: %.2f" % train_error) #Error de test y_predict = model_elasticCV.predict(X_test) elastic_error = mean_squared_error(y_test, y_predict, squared = False) print("RMSE test: %.2f" % elastic_error) RMSE entrenamiento: 4.87 RMSE test: 3.86 El error es parecido a Lasso. ===== Comparaciones entre modelos ===== df_comparacion = pd.DataFrame({ 'modelo': ['OLS', 'Ridge', 'Lasso', 'Lasso_1std', 'Elastic-net'], 'test rmse': [ols_error, ridge_error, lasso_error, lasso_1std_error, elastic_error] }) fig, ax = plt.subplots(figsize=(7, 3.84)) df_comparacion.set_index('modelo').plot(kind='barh', ax=ax) ax.set_xlabel('rmse') ax.set_ylabel('modelo') ax.set_title('Comparación de modelos') {{ :clase:ia:saa:6_optimizacion_modelo:comparacion.png?600 |}} Para este caso, el modelo con el menor error de test es Ridge. Lasso_1std empeora un poco el error en comparación con el modelo simple, pero utilizando sólo 17 predictores. ===== Referencias ===== [[https://www.cienciadedatos.net/documentos/py14-ridge-lasso-elastic-net-python.html|Regularización Ridge, Lasso y Elastic Net con Python by Joaquín Amat Rodrigo, available under a Attribution 4.0 International (CC BY 4.0)]]