Morten Hjorth-Jensen
101 lines
4.5 KiB
Python
101 lines
4.5 KiB
Python
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
|
|
class LassoGD:
|
|
def __init__(self, lr=0.01, l1_penalty=1.0, tol=1e-6, max_iter=1000):
|
|
"""Initialize LASSO regressor with given hyperparameters."""
|
|
self.lr = lr # learning rate (step size)
|
|
self.l1_penalty = l1_penalty # L1 regularization strength (λ)
|
|
self.tol = tol # convergence tolerance for change in cost
|
|
self.max_iter = max_iter # maximum iterations to run
|
|
self.weights = None # model weights (including intercept as w[0])
|
|
self.feature_means = None # to store feature means for normalization
|
|
self.feature_stds = None # to store feature std devs for normalization
|
|
self.cost_history = None # to record cost at each iteration
|
|
|
|
def fit(self, X, y, normalize=True):
|
|
"""Train the LASSO model on data X (shape m x n) and targets y (length m)."""
|
|
X = np.array(X, dtype=float)
|
|
y = np.array(y, dtype=float)
|
|
m, n = X.shape
|
|
# 1. Feature normalization (zero mean, unit variance)
|
|
if normalize:
|
|
self.feature_means = X.mean(axis=0)
|
|
self.feature_stds = X.std(axis=0)
|
|
self.feature_stds[self.feature_stds == 0] = 1.0 # avoid division by zero
|
|
X = (X - self.feature_means) / self.feature_stds
|
|
else:
|
|
# If not normalizing, set means=0 and stds=1 for consistency
|
|
self.feature_means = np.zeros(n)
|
|
self.feature_stds = np.ones(n)
|
|
# Add bias term (intercept) as an extra column of ones in X
|
|
X_bias = np.hstack([np.ones((m, 1)), X])
|
|
# Initialize weights (n features + 1 intercept) to zero
|
|
self.weights = np.zeros(n + 1)
|
|
self.cost_history = []
|
|
|
|
prev_cost = float('inf')
|
|
# Gradient Descent Loop
|
|
for it in range(self.max_iter):
|
|
# 2. Predictions for current weights
|
|
y_pred = X_bias.dot(self.weights)
|
|
error = y_pred - y
|
|
# 3. Compute cost = MSE + L1 penalty (do not penalize intercept w[0])
|
|
mse_cost = (error ** 2).mean() / 2.0
|
|
l1_cost = self.l1_penalty * np.sum(np.abs(self.weights[1:]))
|
|
cost = mse_cost + l1_cost
|
|
self.cost_history.append(cost)
|
|
# Check convergence: stop if change in cost is below tolerance
|
|
if abs(prev_cost - cost) < self.tol:
|
|
break
|
|
prev_cost = cost
|
|
# 4. Compute gradient of MSE part
|
|
grad_mse = (X_bias.T.dot(error)) / m # gradient of 1/(2m)*RSS is X^T(error)/m
|
|
# 5. Perform gradient descent update with L1 penalty via soft-thresholding
|
|
# Take a gradient step for MSE
|
|
w_temp = self.weights - self.lr * grad_mse
|
|
# Soft-thresholding for L1: shrink weights toward 0 by lr*λ
|
|
thresh = self.lr * self.l1_penalty
|
|
w0 = w_temp[0] # intercept (no regularization)
|
|
w_rest = w_temp[1:]
|
|
# Apply soft threshold to each weight in w_rest
|
|
w_rest_updated = np.sign(w_rest) * np.maximum(np.abs(w_rest) - thresh, 0.0)
|
|
# Update weights (combine intercept and rest)
|
|
self.weights = np.concatenate(([w0], w_rest_updated))
|
|
# End of gradient descent loop
|
|
|
|
def predict(self, X):
|
|
"""Make predictions using the trained model on new data X."""
|
|
X = np.array(X, dtype=float)
|
|
# Normalize using the training mean and std
|
|
X_norm = (X - self.feature_means) / self.feature_stds
|
|
# Add bias term
|
|
X_bias = np.hstack([np.ones((X_norm.shape[0], 1)), X_norm])
|
|
return X_bias.dot(self.weights)
|
|
|
|
# --- Example usage on a synthetic dataset ---
|
|
np.random.seed(0)
|
|
# Create synthetic data: m samples, n features
|
|
m, n = 100, 5
|
|
X = np.random.randn(m, n)
|
|
# True underlying weights for features (some are zero to illustrate feature selection)
|
|
true_w = np.array([0, 0, 5, 0, -3], dtype=float)
|
|
true_intercept = 10.0
|
|
# Generate targets with a linear combination of X and noise
|
|
y = true_intercept + X.dot(true_w) + 0.5 * np.random.randn(m)
|
|
|
|
# Train LASSO regression model
|
|
model = LassoGD(lr=0.05, l1_penalty=0.5, tol=1e-6, max_iter=1000)
|
|
model.fit(X, y)
|
|
print("Learned weights (intercept + coefficients):", model.weights)
|
|
|
|
# Plot the cost function history over iterations
|
|
plt.figure(figsize=(6,4))
|
|
plt.plot(model.cost_history, label="Cost")
|
|
plt.title("Cost Function Value vs Iterations")
|
|
plt.xlabel("Iteration")
|
|
plt.ylabel("Cost (MSE + L1 penalty)")
|
|
plt.legend()
|
|
plt.grid(True)
|
|
plt.show()
|