lars/FYS-STK4155
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

update

Browse Source
This commit is contained in:
Morten Hjorth-Jensen
2025-09-01 08:19:41 +02:00
parent cfa272a000
commit 2015cf2c36
3 changed files with 528 additions and 244 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/src/week36/exercisesweek37.do.txt
-167
View File
@@ -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.
doc/src/week37/exercisesweek37.do.txt
+168 -77
View File
@@ -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.