From 2015cf2c36ff63394b1a0dec280f5ccc79bbefc3 Mon Sep 17 00:00:00 2001 From: Morten Hjorth-Jensen Date: Mon, 1 Sep 2025 08:19:41 +0200 Subject: [PATCH] update --- doc/LectureNotes/exercisesweek37.ipynb | 360 +++++++++++++++++++++++++ doc/src/week36/exercisesweek37.do.txt | 167 ------------ doc/src/week37/exercisesweek37.do.txt | 245 +++++++++++------ 3 files changed, 528 insertions(+), 244 deletions(-) create mode 100644 doc/LectureNotes/exercisesweek37.ipynb delete mode 100644 doc/src/week36/exercisesweek37.do.txt diff --git a/doc/LectureNotes/exercisesweek37.ipynb b/doc/LectureNotes/exercisesweek37.ipynb new file mode 100644 index 000000000..68dd47235 --- /dev/null +++ b/doc/LectureNotes/exercisesweek37.ipynb @@ -0,0 +1,360 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "7d56b2d5", + "metadata": { + "editable": true + }, + "source": [ + "\n", + "" + ] + }, + { + "cell_type": "markdown", + "id": "c7a8e9c7", + "metadata": { + "editable": true + }, + "source": [ + "# Exercises week 36\n", + "**Implementing gradient descent for Ridge and ordinary Least Squares Regression**\n", + "\n", + "Date: **September 8-12, 2025**" + ] + }, + { + "cell_type": "markdown", + "id": "cf8f0ecb", + "metadata": { + "editable": true + }, + "source": [ + "## Learning goals\n", + "\n", + "After having completed these exercises you will have:\n", + "1. Your own code for the implementation of the simplest gradient descent approach applied to ordinary least squares (OLS) and Ridge regression\n", + "\n", + "2. Be able to compare the analytical expressions for OLS and Rudge regression with the gradient descent approach\n", + "\n", + "3. Explore the role of the learning rate in the gradient descent approach and the hyperparameter $\\lambda$ in Ridge regression\n", + "\n", + "4. Scale the data properly" + ] + }, + { + "cell_type": "markdown", + "id": "a67ae548", + "metadata": { + "editable": true + }, + "source": [ + "## Ridge regression and a new Synthetic Dataset\n", + "\n", + "We create a synthetic linear regression dataset with a sparse\n", + "underlying relationship. This means we have many features but only a\n", + "few of them actually contribute to the target. In our example, we’ll\n", + "use 10 features with only 3 non-zero weights in the true model. This\n", + "way, the target is generated as a linear combination of a few features\n", + "(with known coefficients) plus some random noise. The steps we include are:\n", + "\n", + "Decide on the number of samples and features (e.g. 100 samples, 10 features).\n", + "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.\n", + "\n", + "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.\n", + "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).\n", + "\n", + "Below is the code to generate the dataset:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "f2d4a55d", + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "# Set random seed for reproducibility\n", + "np.random.seed(0)\n", + "\n", + "# Define dataset size\n", + "n_samples = 100\n", + "n_features = 10\n", + "\n", + "# Define true coefficients (sparse linear relationship)\n", + "theta_true = np.array([5.0, -3.0, 0.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0])\n", + "\n", + "# Generate feature matrix X (n_samples x n_features) with random values\n", + "X = np.random.randn(n_samples, n_features) # standard normal distribution\n", + "\n", + "# Generate target values y with a linear combination of X and theta_true, plus noise\n", + "noise = 0.5 * np.random.randn(n_samples) # Gaussian noise\n", + "y = X.dot @ theta_true + noise" + ] + }, + { + "cell_type": "markdown", + "id": "a445583b", + "metadata": { + "editable": true + }, + "source": [ + "This code produces a dataset where only features 0, 1, and 6\n", + "significantly influence $\\boldsymbol{y}$. The rest of the features have zero true\n", + "coefficient, so they only contribute noise. For example, feature 0 has\n", + "a true weight of 5.0, feature 1 has -3.0, and feature 6 has 2.0, so\n", + "the expected relationship is:" + ] + }, + { + "cell_type": "markdown", + "id": "4a81ddf9", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "y \\approx 5 \\times X_0 \\;-\\; 3 \\times X_1 \\;+\\; 2 \\times X_6 \\;+\\; \\text{noise}.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "ae590275", + "metadata": { + "editable": true + }, + "source": [ + "## Exercise 1, scale your data\n", + "\n", + "Before fitting a regression model, it is good practice to normalize or\n", + "standardize the features. This ensures all features are on a\n", + "comparable scale, which is especially important when using\n", + "regularization. Here we will perform standardization, scaling each\n", + "feature to have mean 0 and standard deviation 1:\n", + "\n", + "Compute the mean and standard deviation of each column (feature) in $bm{X}$.\n", + "Subtract the mean and divide by the standard deviation for each feature.\n", + "\n", + "We will also center the target $\\boldsymbol{y}$ to mean $0$. Centering $\\boldsymbol{y}$\n", + "(and each feature) means the model won’t require a separate intercept\n", + "term – the data is shifted such that the intercept is effectively 0\n", + ". (In practice, one could include an intercept in the model and not\n", + "penalize it, but here we simplify by centering.)" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "8b40c47a", + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "# Standardize features (zero mean, unit variance for each feature)\n", + "X_mean = X.mean(axis=0)\n", + "X_std = X.std(axis=0)\n", + "X_std[X_std == 0] = 1 # safeguard to avoid division by zero for constant features\n", + "X_norm = (X - X_mean) / X_std\n", + "\n", + "# Center the target to zero mean (optional, to simplify intercept handling)\n", + "y_mean = ?\n", + "y_centered = ?" + ] + }, + { + "cell_type": "markdown", + "id": "ff9c0c81", + "metadata": { + "editable": true + }, + "source": [ + "### 1a)\n", + "\n", + "Fill in the necessary details.\n", + "\n", + "After this preprocessing, each column of $\\boldsymbol{X}_norm$ has mean zero and standard deviation $1$\n", + "and $\\boldsymbol{y}_centered$ has mean 0. This makes the optimization landscape\n", + "nicer and ensures the regularization penalty $\\lambda \\sum_j\n", + "\\beta_j^2$ treats each coefficient fairly (since features are on the\n", + "same scale)." + ] + }, + { + "cell_type": "markdown", + "id": "d27c70e4", + "metadata": { + "editable": true + }, + "source": [ + "## Exercise 2, use the analytical formulae for OLS and Ridge regression to find the optimal paramters $\\boldsymbol{theta}$" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "9f1e5184", + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "# Set regularization parameter, either a single value or a vector of values\n", + "lambda = ?\n", + "\n", + "# 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\n", + "I = np.eye(n_features)\n", + "theta_closed_formRidge = ?\n", + "theta_closed_formOLS = ?\n", + "\n", + "print(\"Closed-form Ridge coefficients:\", theta_closed_form)\n", + "print(\"Closed-form OLS coefficients:\", theta_closed_form)" + ] + }, + { + "cell_type": "markdown", + "id": "2ec556b9", + "metadata": { + "editable": true + }, + "source": [ + "This computes the ridge and OLS regression coefficients directly. The identity\n", + "matrix $I$ has the same size as $X^T X$ (which is n_features x\n", + "n_features), and lam * I adds $\\lambda$ to the diagonal of $X^T X. We\n", + "then invert this matrix and multiply by $X^T y. The result\n", + "for $\\boldsymbol{\\theta}$ is a NumPy array of shape (n_features,) containing the\n", + "fitted weights." + ] + }, + { + "cell_type": "markdown", + "id": "a821f0c5", + "metadata": { + "editable": true + }, + "source": [ + "### 2a)\n", + "\n", + "Finalize the OLS and Ridge regression determination of the optimal parameters $bm{\\theta}$." + ] + }, + { + "cell_type": "markdown", + "id": "d637130e", + "metadata": { + "editable": true + }, + "source": [ + "### 2b)\n", + "\n", + "Explore the results as function of different values of the hyperparameter $\\lambda$. See for example exercise 4 from week 36." + ] + }, + { + "cell_type": "markdown", + "id": "b455ce7e", + "metadata": { + "editable": true + }, + "source": [ + "## Implementing the simplest form for gradient descent\n", + "\n", + "Alternatively, we can fit the ridge regression model using gradient\n", + "descent. This is useful to visualize the iterative convergence and is\n", + "necessary if $n$ and $p$ are so large that the closed-form might be\n", + "too slow or memory-intensive. We derive the gradients from the cost\n", + "functions defined above. Use the gradients of the Ridge and OLS cost functions with respect to\n", + "the parameters $\\boldsymbol{\\theta}$ and set up (using the template below) your own gradient descent code for OLS and Ridge regression.\n", + "\n", + "Below is a template code for gradient descent implementation of ridge:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "cfa1eb29", + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "# Gradient descent parameters, learning rate eta first\n", + "eta = 0.1\n", + "# Then number of iterations\n", + "num_iters = 1000\n", + "\n", + "# Initialize weights for gradient descent\n", + "theta = np.zeros(n_features)\n", + "\n", + "# Arrays to store history for plotting\n", + "cost_history = np.zeros(num_iters)\n", + "\n", + "# Gradient descent loop\n", + "m = n_samples # number of examples\n", + "for t in range(num_iters):\n", + " # Compute prediction error\n", + " error = X_norm.dot(theta) - y_centered \n", + " # Compute cost for OLS and Ridge (MSE + regularization for Ridge) for monitoring\n", + " cost_OLS = ?\n", + " cost_Ridge = ?\n", + " cost_history[t] = ?\n", + " # Compute gradients for OSL and Ridge\n", + " grad_OLS = ?\n", + " grad_Ridge = ?\n", + " # Update parameters theta\n", + " theta_gdOLS = ?\n", + " theta_gdRidge = ? \n", + "\n", + "# After the loop, theta contains the fitted coefficients\n", + "theta_gdOLS = ?\n", + "theta_gdRidge = ?\n", + "print(\"Gradient Descent OLS coefficients:\", theta_gdOLS)\n", + "print(\"Gradient Descent Ridge coefficients:\", theta_gdRidge)" + ] + }, + { + "cell_type": "markdown", + "id": "dc78d58d", + "metadata": { + "editable": true + }, + "source": [ + "### 3a)\n", + "\n", + "Discuss the results as function of the learning rate paramaters and the number of iterations." + ] + }, + { + "cell_type": "markdown", + "id": "15060acb", + "metadata": { + "editable": true + }, + "source": [ + "### 3b)\n", + "\n", + "Add a stopping parameter as function of the number iterations. \n", + "\n", + "If everything worked correctly, the learned coefficients should be\n", + "close to the true values [5.0, -3.0, 0.0, …, 2.0, …] that we used to\n", + "generate the data. Keep in mind that due to regularization and noise,\n", + "the learned values will not exactly equal the true ones, but they\n", + "should be in the same ballpark." + ] + } + ], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/doc/src/week36/exercisesweek37.do.txt b/doc/src/week36/exercisesweek37.do.txt deleted file mode 100644 index 114fec64d..000000000 --- a/doc/src/week36/exercisesweek37.do.txt +++ /dev/null @@ -1,167 +0,0 @@ -TITLE: Week 37: Linear Regression and Gradient descent -AUTHOR: Morten Hjorth-Jensen {copyright, 1999-present|CC BY-NC} at Department of Physics, University of Oslo, Norway -DATE: September 1-5, 2025 - - - -===== 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{\bm{\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 $\bm{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 $\bm{X}\hat{\bm{\theta}}$ and add some noise (to simulate measurement error or unexplained variance). - - -Below is the code to generate the dataset: -!bc pycod -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 -!ec - -This code produces a dataset where only features 0, 1, and 6 -significantly influence $\bm{y}$. The rest of the features have zero true -coefficient, so they only contribute noise. 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: -!bt -\[ -y \approx 5 \times X_0 \;-\; 3 \times X_1 \;+\; 2 \times X_6 \;+\; \text{noise}. -\] -!et - - - -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: - -Compute the mean and standard deviation of each column (feature) in $bm{X}$. -Subtract the mean and divide by the standard deviation for each feature. - - -We also center the target $\bm{y}$ to mean $0$. Centering $\bm{y}$ -(and each feature) means the model won’t 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.) - -!bc pycod -# 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() -y_centered = y - y_mean -!ec - -After this preprocessing, each column of $\bm{X}_norm$ has mean zero and standard deviation $1$ -and $\bm{y}_centered$ has mean 0. This makes the optimization landscape -nicer and ensures the regularization penalty $\lambda \sum_j -\beta_j^2$ treats each coefficient fairly (since features are on the -same scale). - -!bc pycod -# Set regularization parameter -lam = 1.0 - -# Closed-form Ridge solution: w = (X^T X + lam * I)^{-1} X^T y -I = np.eye(n_features) -w_closed_form = np.linalg.inv(X_norm.T.dot(X_norm) + lam * I).dot(X_norm.T).dot(y_centered) - -print("Closed-form Ridge coefficients:", w_closed_form) -!ec - -This computes the ridge regression coefficients directly. The identity -matrix $I$ has the same size as $X^T X$ (which is n_features x -n_features), and lam * I adds $\lambda$ to the diagonal of $X^T X. We -then invert this matrix and multiply by $X^T y. The result -for $\bm{\theta}$ is a NumPy array of shape (n_features,) containing the -fitted weights. - - - -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 -function defined above. The gradient of the ridge cost with respect to -the weight vector $w$ is: - - - -Below is the code for gradient descent implementation of ridge: -!bc pycod -# Gradient descent parameters -alpha = 0.1 -num_iters = 1000 - -# Initialize weights for gradient descent -theta = np.zeros(n_features) - -# Arrays to store history for plotting -cost_history = np.zeros(num_iters) - -# Gradient descent loop -m = n_samples # number of examples -for t in range(num_iters): - # Compute prediction error - error = X_norm.dot(theta) - y_centered # shape (m,) - # Compute cost (MSE + regularization) for monitoring - cost = (1/(2*m)) * np.dot(error, error) + (lam/(2*m)) * np.dot(theta, theta) - cost_history[t] = cost - # Compute gradient - grad = (1/m) * (X_norm.T.dot(error) + lam * theta) - # Update weights - theta = theta - alpha * grad - -# After the loop, theta contains the fitted coefficients -theta_gd = theta -print("Gradient Descent Ridge coefficients:", theta_gd) -!ec - - -Let us confirm that the two approaches (closed-form and gradient -descent) give similar results, and then evaluate the model. First, -compare the learned coefficients to the true coefficients: - -!bc pycod -print("True coefficients:", theta_true) -print("Closed-form learned coefficients:", theta_closed_form) -print("Gradient descent learned coefficients:", theta_gd) -!ec -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. - diff --git a/doc/src/week37/exercisesweek37.do.txt b/doc/src/week37/exercisesweek37.do.txt index 9950aa9d0..0a19d1538 100644 --- a/doc/src/week37/exercisesweek37.do.txt +++ b/doc/src/week37/exercisesweek37.do.txt @@ -1,102 +1,193 @@ -TITLE: Exercises week 37 -AUTHOR: September 9-13, 2024 -DATE: Deadline is Friday September 13 at midnight +TITLE: Exercises week 36 +AUTHOR: Implementing gradient descent for Ridge and ordinary Least Squares Regression +DATE: September 8-12, 2025 -===== Overarching aims of the exercises this week ===== +===== Learning goals ===== + +After having completed these exercises you will have: +o Your own code for the implementation of the simplest gradient descent approach applied to ordinary least squares (OLS) and Ridge regression +o Be able to compare the analytical expressions for OLS and Rudge regression with the gradient descent approach +o Explore the role of the learning rate in the gradient descent approach and the hyperparameter $\lambda$ in Ridge regression +o Scale the data properly + +===== Ridge regression and a new Synthetic Dataset ===== -This exercise deals with various mean values and variances in linear -regression method (here it may be useful to look up chapter 3, -equation (3.8) of "Trevor Hastie, Robert Tibshirani, Jerome -H. Friedman, The Elements of Statistical Learning, -Springer":"https://www.springer.com/gp/book/9780387848570"). The -exercise is also a part of project 1 and can be reused in the theory -part of the project. +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: -For more discussions on Ridge regression and calculation of -expectation values, "Wessel van -Wieringen's":"https://arxiv.org/abs/1509.09169" article is highly -recommended. +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{\bm{\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 $\bm{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 $\bm{X}\hat{\bm{\theta}}$ and add some noise (to simulate measurement error or unexplained variance). -The assumption we have made is that there exists a continuous function -$f(\bm{x})$ and a normal distributed error $\bm{\varepsilon}\sim N(0, -\sigma^2)$ which describes our data +Below is the code to generate the dataset: +!bc pycod +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 +!ec + +This code produces a dataset where only features 0, 1, and 6 +significantly influence $\bm{y}$. The rest of the features have zero true +coefficient, so they only contribute noise. 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: !bt \[ -\bm{y} = f(\bm{x})+\bm{\varepsilon} -\] -!et - -We then approximate this function $f(\bm{x})$ with our model $\bm{\tilde{y}}$ from the solution of the linear regression equations (ordinary least squares OLS), that is our -function $f$ is approximated by $\bm{\tilde{y}}$ where we minimized $(\bm{y}-\bm{\tilde{y}})^2$, with -!bt -\[ -\bm{\tilde{y}} = \bm{X}\bm{\beta}. -\] -!et -The matrix $\bm{X}$ is the so-called design or feature matrix. - -===== Exercise: Expectation values for ordinary least squares expressions ===== - -Show that the expectation value of $\bm{y}$ for a given element $i$ -!bt -\[ -\mathbb{E}(y_i) =\sum_{j}x_{ij} \beta_j=\mathbf{X}_{i, \ast} \, \bm{\beta}, -\] -!et -and that -its variance is -!bt -\[ -\mbox{Var}(y_i) = \sigma^2. -\] -!et -Hence, $y_i \sim N( \mathbf{X}_{i, \ast} \, \bm{\beta}, \sigma^2)$, that is $\bm{y}$ follows a normal distribution with -mean value $\bm{X}\bm{\beta}$ and variance $\sigma^2$. - -With the OLS expressions for the optimal parameters $\bm{\hat{\beta}}$ show that -!bt -\[ -\mathbb{E}(\bm{\hat{\beta}}) = \bm{\beta}. -\] -!et -Show finally that the variance of $\bm{\bm{\beta}}$ is -!bt -\[ -\mbox{Var}(\bm{\hat{\beta}}) = \sigma^2 \, (\mathbf{X}^{T} \mathbf{X})^{-1}. +y \approx 5 \times X_0 \;-\; 3 \times X_1 \;+\; 2 \times X_6 \;+\; \text{noise}. \] !et -We can use the last expression when we define a "so-called confidence interval":"https://en.wikipedia.org/wiki/Confidence_interval" for the parameters $\beta$. -A given parameter $\beta_j$ is given by the diagonal matrix element of the above matrix. +===== 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: + +Compute the mean and standard deviation of each column (feature) in $bm{X}$. +Subtract the mean and divide by the standard deviation for each feature. -===== Exercise: Expectation values for Ridge regression ===== +We will also center the target $\bm{y}$ to mean $0$. Centering $\bm{y}$ +(and each feature) means the model won’t 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.) -Show that -!bt -\[ -\mathbb{E} \big[ \hat{\bm{\beta}}^{\mathrm{Ridge}} \big]=(\mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I}_{pp})^{-1} (\mathbf{X}^{\top} \mathbf{X})\bm{\beta}. -\] -!et -We see clearly that -$\mathbb{E} \big[ \hat{\bm{\beta}}^{\mathrm{Ridge}} \big] \not= \mathbb{E} \big[\hat{\bm{\beta}}^{\mathrm{OLS}}\big ]$ for any $\lambda > 0$. +!bc pycod +# 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_centered = ? +!ec + +=== 1a) === +Fill in the necessary details. + +After this preprocessing, each column of $\bm{X}_norm$ has mean zero and standard deviation $1$ +and $\bm{y}_centered$ has mean 0. This makes the optimization landscape +nicer and ensures the regularization penalty $\lambda \sum_j +\beta_j^2$ treats each coefficient fairly (since features are on the +same scale). -Show also that the variance is +===== Exercise 2, use the analytical formulae for OLS and Ridge regression to find the optimal paramters $\bm{theta}$ ===== -!bt -\[ -\mbox{Var}[\hat{\bm{\beta}}^{\mathrm{Ridge}}]=\sigma^2[ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1} \mathbf{X}^{T}\mathbf{X} \{ [ \mathbf{X}^{\top} \mathbf{X} + \lambda \mathbf{I} ]^{-1}\}^{T}, -\] -!et -and it is easy to see that if the parameter $\lambda$ goes to infinity then the variance of the Ridge parameters $\bm{\beta}$ goes to zero. +!bc pycod +# Set regularization parameter, either a single value or a vector of values +lambda = ? + +# 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 = ? +theta_closed_formOLS = ? + +print("Closed-form Ridge coefficients:", theta_closed_form) +print("Closed-form OLS coefficients:", theta_closed_form) +!ec + +This computes the ridge and OLS regression coefficients directly. The identity +matrix $I$ has the same size as $X^T X$ (which is n_features x +n_features), and lam * I adds $\lambda$ to the diagonal of $X^T X. We +then invert this matrix and multiply by $X^T y. The result +for $\bm{\theta}$ is a NumPy array of shape (n_features,) containing the +fitted weights. + +=== 2a) === +Finalize the OLS and Ridge regression determination of the optimal parameters $bm{\theta}$. + +=== 2b) === +Explore the results as function of different values of the hyperparameter $\lambda$. See for example exercise 4 from week 36. + +===== 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 $\bm{\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: +!bc pycod +# Gradient descent parameters, learning rate eta first +eta = 0.1 +# Then number of iterations +num_iters = 1000 + +# Initialize weights for gradient descent +theta = np.zeros(n_features) + +# Arrays to store history for plotting +cost_history = np.zeros(num_iters) + +# Gradient descent loop +m = n_samples # number of examples +for t in range(num_iters): + # Compute prediction error + error = X_norm.dot(theta) - y_centered + # Compute cost for OLS and Ridge (MSE + regularization for Ridge) for monitoring + cost_OLS = ? + cost_Ridge = ? + cost_history[t] = ? + # Compute gradients for OSL and Ridge + grad_OLS = ? + grad_Ridge = ? + # Update parameters theta + theta_gdOLS = ? + theta_gdRidge = ? + +# After the loop, theta contains the fitted coefficients +theta_gdOLS = ? +theta_gdRidge = ? +print("Gradient Descent OLS coefficients:", theta_gdOLS) +print("Gradient Descent Ridge coefficients:", theta_gdRidge) +!ec + +=== 3a) === +Discuss the results as function of the learning rate paramaters and the number of iterations. + +=== 3b) === +Add a stopping parameter as function of the number iterations. + +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. +