Exercises week 37¶

Implementing gradient descent for Ridge and ordinary Least Squares Regression

Date: September 8-12, 2025

Python Code can be found at https://github.uio.no/larsbog/FYS-STK4155

Learning goals¶

After having completed these exercises you will have:

  1. Your own code for the implementation of the simplest gradient descent approach applied to ordinary least squares (OLS) and Ridge regression

  2. Be able to compare the analytical expressions for OLS and Rudge regression with the gradient descent approach

  3. Explore the role of the learning rate in the gradient descent approach and the hyperparameter $\lambda$ in Ridge regression

  4. Scale the data properly

Simple one-dimensional second-order polynomial¶

We start with a very simple function

$$ f(x)= 2-x+5x^2, $$

defined for $x\in [-2,2]$. You can add noise if you wish.

We are going to fit this function with a polynomial ansatz. The easiest thing is to set up a second-order polynomial and see if you can fit the above function. Feel free to play around with higher-order polynomials.

Exercise 1, scale your data¶

Before fitting a regression model, it is good practice to normalize or standardize the features. This ensures all features are on a comparable scale, which is especially important when using regularization. Here we will perform standardization, scaling each feature to have mean 0 and standard deviation 1.

1a)¶

Compute the mean and standard deviation of each column (feature) in your design/feature matrix $\boldsymbol{X}$. Subtract the mean and divide by the standard deviation for each feature.

We will also center the target $\boldsymbol{y}$ to mean $0$. Centering $\boldsymbol{y}$ (and each feature) means the model does not require a separate intercept term, the data is shifted such that the intercept is effectively 0 . (In practice, one could include an intercept in the model and not penalize it, but here we simplify by centering.) Choose $n=100$ data points and set up $\boldsymbol{x}, $\boldsymbol{y}$ and the design matrix $\boldsymbol{X}$.

In [1]:
import numpy as np


def polynomial_features(x, p, intercept=False):
    n = len(x)
    if intercept:
        P = np.arange(p+1)
    else:
        P = np.arange(1, p+1)
    X = np.power(x[:, np.newaxis], P)
    return X


theta_true = np.array([-2.0, 5.0])
X = polynomial_features(np.linspace(0, 1, 100), p=2)
y = 2 + theta_true@X.T + 0.1 * np.random.randn(100)
In [2]:
# Standardize features (zero mean, unit variance for each feature)
X_mean = X.mean(axis=0)
X_std = X.std(axis=0)
X_std[X_std == 0] = 1  # safeguard to avoid division by zero for constant features
X_norm = (X - X_mean) / X_std

# Center the target to zero mean (optional, to simplify intercept handling)
y_mean = y.mean(axis=0)
y_centered = y - y_mean

n_features = X_norm.shape[1]
theta_true *= X_std
In [3]:
theta_true
Out[3]:
array([-0.58315293,  1.50663026])
In [4]:
X_std
Out[4]:
array([0.29157647, 0.30132605])
In [5]:
import matplotlib.pyplot as plt

plt.scatter(X_norm[:, 0], y_centered)
plt.xlabel("Feature 1")
plt.ylabel("Target")
Out[5]:
Text(0, 0.5, 'Target')
No description has been provided for this image

Fill in the necessary details.

After this preprocessing, each column of $\boldsymbol{X}_{\mathrm{norm}}$ has mean zero and standard deviation $1$ and $\boldsymbol{y}_{\mathrm{centered}}$ has mean 0. This makes the optimization landscape nicer and ensures the regularization penalty $\lambda \sum_j \theta_j^2$ in Ridge regression treats each coefficient fairly (since features are on the same scale).

Exercise 2, calculate the gradients¶

Find the gradients for OLS and Ridge regression using the mean-squared error as cost/loss function.

The gradients can be calculated as

$$ \nabla_\theta C_{OLS} = \nabla_\theta (X\cdot\theta - y)^2 = 2 X (X\cdot \theta - y) $$

The gradient of the Ridge cost function then directly follows as

$$ \nabla_\theta C_\mathrm{Ridge} = \nabla_\theta C_{OLS} + \nabla_\theta \lambda \theta^2 = 2 X (X\cdot \theta - y) + 2 \lambda \theta $$

Because the factors of 2 are somewhat tedious, I will from here on use $\tilde C = \frac{C}{2}$.

Exercise 3, use the analytical formulae for OLS and Ridge regression to find the optimal paramters $\boldsymbol{\theta}$¶

In [6]:
def Ridge_parameters(X, y, lam = 0.01):
    # Assumes X is scaled and has no intercept column
    return np.linalg.inv(X.T @ X + lam * np.eye(X.shape[1])) @ X.T @ y
def OLS_parameters(X, y):
    return Ridge_parameters(X, y, lam = 0.0)
In [7]:
# Set regularization parameter, either a single value or a vector of values
lambda_ = 0.1

# Analytical form for OLS and Ridge solution: theta_Ridge = (X^T X + lambda * I)^{-1} X^T y and theta_OLS = (X^T X)^{-1} X^T y
I = np.eye(n_features)
theta_closed_formRidge = Ridge_parameters(X_norm, y_centered, lam=lambda_)
theta_closed_formOLS = OLS_parameters(X_norm, y_centered)

print("Closed-form Ridge coefficients:", theta_closed_formRidge)
print("Closed-form OLS coefficients:", theta_closed_formOLS)

Closed-form Ridge coefficients: [-0.4835408   1.41698898]
Closed-form OLS coefficients: [-0.51267301  1.44659558]

This computes the Ridge and OLS regression coefficients directly. The identity matrix $I$ has the same size as $X^T X$. It adds $\lambda$ to the diagonal of $X^T X$ for Ridge regression. We then invert this matrix and multiply by $X^T y$. The result for $\boldsymbol{\theta}$ is a NumPy array of shape (n$\_$features,) containing the fitted parameters $\boldsymbol{\theta}$.

3a)¶

Finalize, in the above code, the OLS and Ridge regression determination of the optimal parameters $\boldsymbol{\theta}$.

3b)¶

Explore the results as function of different values of the hyperparameter $\lambda$. See for example exercise 4 from week 36.

In [8]:
import matplotlib.pyplot as plt

n_lam = 8
thetas = np.zeros((n_features, n_lam))
lambdas = np.concatenate([[0], np.logspace(-5, 0, n_lam-1)])


for i, lam in enumerate(lambdas):
    thetas[:, i] = Ridge_parameters(X_norm, y_centered, lam=lam)

avg_thetas = np.mean(thetas, axis=1)
norm_thetas = thetas - avg_thetas[:, np.newaxis]

fig, ax = plt.subplots(figsize=(8,6))
# Annotation
im = ax.imshow(norm_thetas, aspect='auto', cmap='viridis')
for i in range(n_features):
    for j in range(n_lam):
        text = ax.text(j, i, f"{norm_thetas[i, j]:.2e}",
                       ha="center", va="center", color="w" if abs(norm_thetas[i, j]) < 0.5 else "black")
ax.set_yticks(np.arange(n_features))
ax.set_xticks(np.arange(n_lam))
ax.set_xticklabels([f"{l:.2e}" for l in lambdas], rotation=45)
ax.set_ylabel("Features")
ax.set_xlabel("Regularization Parameter (lambda)")
ax.set_title(r"Feature Coefficients ($\theta^\lambda_i - \bar{\theta}_i$)")
Out[8]:
Text(0.5, 1.0, 'Feature Coefficients ($\\theta^\\lambda_i - \\bar{\\theta}_i$)')
No description has been provided for this image

Exercise 4, Implementing the simplest form for gradient descent¶

Alternatively, we can fit the ridge regression model using gradient descent. This is useful to visualize the iterative convergence and is necessary if $n$ and $p$ are so large that the closed-form might be too slow or memory-intensive. We derive the gradients from the cost functions defined above. Use the gradients of the Ridge and OLS cost functions with respect to the parameters $\boldsymbol{\theta}$ and set up (using the template below) your own gradient descent code for OLS and Ridge regression.

Below is a template code for gradient descent implementation of ridge:

In [9]:
def gradient_descent(X, y, cost_func, grad_cost_func, eta=0.1, num_iters=1000, **kwargs):
    # Initialize weights
    theta = np.zeros(X.shape[1])
    # Store cost history
    cost_history = np.zeros(num_iters)
    for t in range(num_iters):
        # Compute cost
        cost_history[t] = cost_func(X, y, theta, **kwargs)
        # Compute gradient
        grad = grad_cost_func(X, y, theta, **kwargs)
        # Update weights
        theta -= eta * grad
    return theta, cost_history

def OLS_cost_func(X, y, theta):
    error = X.dot(theta) - y
    return 0.5 * np.mean(error**2)

def OLS_grad_cost_func(X, y, theta):
    error = X.dot(theta) - y
    return X.T.dot(error) / len(y)

def Ridge_cost_func(X, y, theta, lam):
    return OLS_cost_func(X, y, theta) + 0.5 * lam * np.sum(theta**2)

def Ridge_grad_cost_func(X, y, theta, lam):
    return OLS_grad_cost_func(X, y, theta) + lam * theta
In [10]:
eta = 0.1
num_iters = 1000
lam = 1e-3

theta_gdOLS, history_gdOLS = gradient_descent(X_norm, y_centered, OLS_cost_func, OLS_grad_cost_func, eta=eta, num_iters=num_iters)
theta_gdRidge, history_gdRidge = gradient_descent(X_norm, y_centered, Ridge_cost_func, Ridge_grad_cost_func, eta=eta, num_iters=num_iters, lam=lam)

print("Gradient Descent OLS coefficients:", theta_gdOLS)
print("Gradient Descent Ridge coefficients:", theta_gdRidge)

Gradient Descent OLS coefficients: [-0.47433764  1.40826022]
Gradient Descent Ridge coefficients: [-0.44990454  1.38335272]
In [11]:
def theta_error(theta_est):
    return theta_est - theta_true, np.mean((theta_est - theta_true)**2)

print(f"Gradient Descent OLS error: {theta_error(theta_gdOLS)}")
print(f"Gradient Descent Ridge error: {theta_error(theta_gdRidge)}")

Gradient Descent OLS error: (array([ 0.10881529, -0.09837005]), np.float64(0.010758716397832158))
Gradient Descent Ridge error: (array([ 0.13324839, -0.12327754]), np.float64(0.016476242667611243))

4a)¶

Discuss the results as function of the learning rate parameters and the number of iterations.

In [12]:
etas = np.logspace(-4, 0, 10)
max_iterations = np.logspace(1, 5, 10, dtype=int)
theta_errors = np.zeros((len(etas), len(max_iterations), 2))

for i, eta in enumerate(etas):
    for j, max_iter in enumerate(max_iterations):
        # print(f"Running GD with eta={eta}, max_iter={max_iter}")
        theta_gdOLS, _ = gradient_descent(X_norm, y_centered, OLS_cost_func, OLS_grad_cost_func, eta=eta, num_iters=max_iter)
        theta_gdRidge, _ = gradient_descent(X_norm, y_centered, Ridge_cost_func, Ridge_grad_cost_func, eta=eta, num_iters=max_iter, lam=lam)

        theta_errors[i, j, 0] = theta_error(theta_gdOLS)[1]
        theta_errors[i, j, 1] = theta_error(theta_gdRidge)[1]
In [13]:
import matplotlib.cm as mcm
import matplotlib.colors as mcolors

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# Log Colormap
cm = plt.get_cmap('viridis')
# Create LogNorm between min and max of theta_errors
norm = mcolors.LogNorm(vmin=theta_errors.min(), vmax=theta_errors.max())


im1 = ax1.imshow(theta_errors[:, :, 0], aspect='auto', interpolation='nearest', cmap=cm, norm=norm)
ax1.set_title('Gradient Descent OLS Error')
fig.colorbar(im1, ax=ax1)

im2 = ax2.imshow(theta_errors[:, :, 1], aspect='auto', interpolation='nearest', cmap=cm, norm=norm)
ax2.set_title('Gradient Descent Ridge Error')
fig.colorbar(im2, ax=ax2)

for ax in (ax1, ax2):
    ax.set_xlabel('Max Iterations')
    ax.set_ylabel('Learning Rate')
    ax.set_xticks(np.arange(len(max_iterations)))
    ax.set_yticks(np.arange(len(etas)))
    ax.set_xticklabels(max_iterations, rotation=45)
    ax.set_yticklabels(np.round(etas, 4))
No description has been provided for this image
In this simple example (well converging function) the gradient descent converges for all learning rates given enough iterations. The Error on the paramters approaches 10^-3 and stays at this plateau indicating the parameters have converged

4b)¶

Try to add a stopping parameter as function of the number iterations. How would you define a stopping criterion?

We define the stopping criterion via the decrease in the cost function. If the decrease gets too small we stop. This will be implemented using a custom stopping_criterion function so it can easily be replaced.
In [14]:
def gradient_descent(X, y, cost_func, grad_cost_func, eta=0.1, num_iters=1000, stopping_criterion=None, **kwargs):
    # Initialize weights
    theta = np.zeros(X.shape[1])
    # Store cost history
    cost_history = np.zeros(num_iters)
    for t in range(num_iters):
        # Compute cost
        cost_history[t] = cost_func(X, y, theta, **kwargs)
        # Compute gradient
        grad = grad_cost_func(X, y, theta, **kwargs)
        # Update weights
        theta -= eta * grad
        if stopping_criterion is not None and stopping_criterion(cost_history[:t+1]):
            print(f"Converged at iteration {t}")
            break
    return theta, cost_history

def stopping_criterion(cost_history, tol=1e-5):
    if len(cost_history) < 2:
        return False
    return np.abs(cost_history[-1] - cost_history[-2]) < tol

Exercise 5, Ridge regression and a new Synthetic Dataset¶

We create a synthetic linear regression dataset with a sparse underlying relationship. This means we have many features but only a few of them actually contribute to the target. In our example, we’ll use 10 features with only 3 non-zero weights in the true model. This way, the target is generated as a linear combination of a few features (with known coefficients) plus some random noise. The steps we include are:

Decide on the number of samples and features (e.g. 100 samples, 10 features). Define the true coefficient vector with mostly zeros (for sparsity). For example, we set $\hat{\boldsymbol{\theta}} = [5.0, -3.0, 0.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0]$, meaning only features 0, 1, and 6 have a real effect on y.

Then we sample feature values for $\boldsymbol{X}$ randomly (e.g. from a normal distribution). We use a normal distribution so features are roughly centered around 0. Then we compute the target values $y$ using the linear combination $\boldsymbol{X}\hat{\boldsymbol{\theta}}$ and add some noise (to simulate measurement error or unexplained variance).

Below is the code to generate the dataset:

In [15]:
import numpy as np

# Set random seed for reproducibility
np.random.seed(0)

# Define dataset size
n_samples = 100
n_features = 10

# Define true coefficients (sparse linear relationship)
theta_true = np.array([5.0, -3.0, 0.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0])

# Generate feature matrix X (n_samples x n_features) with random values
X = np.random.randn(n_samples, n_features)  # standard normal distribution

# Generate target values y with a linear combination of X and theta_true, plus noise
noise = 0.5 * np.random.randn(n_samples)    # Gaussian noise
y = X.dot(theta_true) + noise

This code produces a dataset where only features 0, 1, and 6 significantly influence $\boldsymbol{y}$. The rest of the features have zero true coefficient. For example, feature 0 has a true weight of 5.0, feature 1 has -3.0, and feature 6 has 2.0, so the expected relationship is:

$$ y \approx 5 \times x_0 \;-\; 3 \times x_1 \;+\; 2 \times x_6 \;+\; \text{noise}. $$

You can remove the noise if you wish to.

Try to fit the above data set using OLS and Ridge regression with the analytical expressions and your own gradient descent codes.

If everything worked correctly, the learned coefficients should be close to the true values [5.0, -3.0, 0.0, …, 2.0, …] that we used to generate the data. Keep in mind that due to regularization and noise, the learned values will not exactly equal the true ones, but they should be in the same ballpark. Which method (OLS or Ridge) gives the best results?

In [16]:
res = gradient_descent(X, y, OLS_cost_func, OLS_grad_cost_func, num_iters=100000, stopping_criterion=stopping_criterion)
print(f"Terminated with error on theta of {theta_error(res[0])[1]:.3e}")

Converged at iteration 85
Terminated with error on theta of 3.750e-03