Ridge-Lasso Regression
In this post, We will review the way of generalization, especially on Ridge and Lasso.
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import lasso_path, Lasso, Ridge
np.random.seed(1)
In logistic Regression, we can define the loss function like this,
$$ L(\beta; \lambda) = -\frac{1}{n} \sum_{i=1}^n (y_i x_i^T \beta - \log(1+ \exp(x_i^T \beta))) $$
For the Regularization, we can add the penalty term. In this case, we added Ridge penalty (also known as L2 error). That is, we measure the squared beta for the penalty.
$$ L(\beta; \lambda) = -\frac{1}{n} \sum_{i=1}^n (y_i x_i^T \beta - \log(1+ \exp(x_i^T \beta))) + \lambda \Vert \beta \Vert_2^2 $$
p = 2
lambda_v = 1.5
n = 10
true_beta = np.array([[1], [-0.5]])
So we need to make predictor. In this example, we assume that the sample data is from binomial distribution.
x = np.random.normal(0, 1, (n, p))
prob = 1 / (1 + np.exp(-x @ true_beta))
prob = prob.reshape((n,))
prob
y = np.random.binomial(np.ones(n, dtype='int32'), prob, n)
y = y.reshape((n, 1))
y
Currently, we generate the ture data from true beta. So how can we get true beta from initial beta?
beta = np.array([.5, .5]).reshape((p, 1))
prob = 1 / (1 + np.exp(-x @ beta))
prob
With Newton raphson method, we differentiate it two times.
$$ g = \frac{\partial}{\partial\beta} L(\beta; \lambda) = -\frac{1}{n} \sum_{i=1}^n \big( y_i x_i - \frac{\exp(x_i^T \beta)}{1 + \exp(x_i^T \beta)} x_i \big) + 2 \lambda \beta \\ H = \frac{\partial^2}{\partial \beta^2} L(\beta ; \lambda) = \frac{1}{n} \sum_{i=1}^n \big( \frac{\exp(x_i^T \beta)}{1 + \exp(x_i^T \beta)} \cdot \frac{1}{1 + \exp(x_i^T \beta)} \big) x_i x_i^T + 2 \lambda I$$grad = np.mean((prob - y) * x, axis=0, keepdims=True).T + 2 * lambda_v * beta
grad
D = np.diag((prob * (1 - prob)).reshape(n))
D[:5, :5]
H = x.T @ D @ x / n + np.diag(np.repeat(2 * lambda_v, p))
H
After we calculate the gradient and hessian matrix of one optimization, we do the one step for the beta update. All we need to do is to repeat these steps until the threshold is reached.
beta = np.zeros((p, 1))
beta_0 = []
beta_1 = []
for i in range(10):
prob = 1 / (1 + np.exp(-x @ beta))
grad = np.mean((prob - y) * x, axis=0, keepdims=True).T + 2 * lambda_v * beta
D = np.diag((prob * (1 - prob)).reshape(n))
H = x.T @ D @ x / n + np.diag(np.repeat(2 * lambda_v, p))
beta_new = beta - np.linalg.inv(H) @ grad
beta_0.append(beta_new[0])
beta_1.append(beta_new[1])
if np.sum(np.abs(beta_new - beta)) < 1e-8:
beta = beta_new
print('Iteration {} beta:'.format(i))
print(beta, '\n')
break
else:
beta = beta_new
print('Iteration {} beta:'.format(i))
print(beta, '\n')
We almost found the solution for beta. So does this solution satisfy the optimality? Or Is the gradient of 0 at this beta? One way to check this is to apply Karush-Kuhn-Tucker (KKT) condition.
prob = 1 / (1 + np.exp(-x @ beta))
grad = np.mean((prob - y) * x, axis=0, keepdims=True).T + 2 * lambda_v * beta
grad
np.all(np.abs(grad) < 1e-8)
In this time, we will use Lasso and Ridge Regression implemented in scikit-learn.
At first, we need to find solution path for each lambda.
Note that, Ridge penalty can be expressed like this,
$$ \beta(\lambda) = \arg\min_{\beta} \frac{1}{2} \Vert Y - X \beta \Vert^2 + \lambda \Vert \beta \Vert^2 $$
n_lambdas = 50
lambda_vec = np.linspace(0, 100, n_lambdas)
coefs = []
for lambda_v in lambda_vec:
ridge = Ridge(alpha=lambda_v, fit_intercept=False)
ridge.fit(x, y)
coefs.append(ridge.coef_)
coefs = np.squeeze(np.array(coefs))
coefs[:6, :]
plt.figure(figsize=(16, 10))
colors = ['b', 'r', 'g']
lstyles = ['-', '--', '-.', ':']
for i in range(p):
l = plt.plot(lambda_vec, coefs[:, i], linestyle=lstyles[i], c=colors[i])
plt.xscale('log')
plt.axis('tight')
plt.xlabel('lambda')
plt.ylabel('beta')
plt.title('Ridge Regression Solution Path')
plt.show()
As you can see, when the $\lambda$ is increased, $\beta$ is closed to 0, which means the 'shrinkage' of the beta.
Unlike Ridge, Lasso uses absolute beta for the penalty.
$$ \beta(\lambda) = \arg\min_{\beta} \frac{1}{2} \Vert Y - X \beta \Vert^2 + \lambda \vert \beta \vert $$
eps = 5e-3
lambdas_lasso, coefs_lasso, _ = lasso_path(x, y, eps=eps, fit_intercept=False)
coefs_lasso = np.squeeze(coefs_lasso)
coefs_lasso[:, :5]
plt.figure(figsize=(16, 10))
colors = ['b', 'r', 'g']
lstyles = ['-', '--', '-.', ':']
for coef, c, l in zip(coefs_lasso, colors, lstyles):
l = plt.plot(lambdas_lasso, coef, linestyle=l, c=c)
plt.axis('tight')
plt.xlabel('lambda')
plt.ylabel('beta')
plt.title('Lasso Regression Solution Path')
plt.show()
