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 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)

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)

El rango de nuestros coeficientes es [-30000, 30000].

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 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)

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)

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)

#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 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)

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)

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)

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).

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 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)

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_)

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')

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)

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.

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) 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)

# 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.

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')

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.

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)