update
This commit is contained in:
Morten Hjorth-Jensen
Binary file not shown.
Binary file not shown.
File diff suppressed because it is too large
Load Diff
File diff suppressed because one or more lines are too long
@@ -419,7 +419,6 @@ document.write(`
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-hessian-matrix">The Hessian matrix</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#simple-program">Simple program</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#id1">Gradient Descent Example</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#and-a-corresponding-example-using-scikit-learn">And a corresponding example using <strong>scikit-learn</strong></a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#gradient-descent-and-ridge">Gradient descent and Ridge</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-hessian-matrix-for-ridge-regression">The Hessian matrix for Ridge Regression</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#program-example-for-gradient-descent-with-ridge-regression">Program example for gradient descent with Ridge Regression</a></li>
|
||||
@@ -631,7 +630,7 @@ This is equivalent to saying that the matrix <span class="math notranslate nohig
|
||||
\]</div>
|
||||
<p>has linearly dependent column vectors, we will not be able to compute the inverse
|
||||
of <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span> and we cannot find the parameters (estimators) <span class="math notranslate nohighlight">\(\theta_i\)</span>.
|
||||
The estimators are only well-defined if <span class="math notranslate nohighlight">\((\boldsymbol{X}^{T}\boldsymbol{X})^{-1}\)</span> exits.
|
||||
The estimators are only well-defined if <span class="math notranslate nohighlight">\((\boldsymbol{X}^{T}\boldsymbol{X})^{-1}\)</span> exists.
|
||||
This is more likely to happen when the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> is high-dimensional. In this case it is likely to encounter a situation where
|
||||
the regression parameters <span class="math notranslate nohighlight">\(\theta_i\)</span> cannot be estimated.</p>
|
||||
<p>A cheap <em>ad hoc</em> approach is simply to add a small diagonal component to the matrix to invert, that is we change</p>
|
||||
@@ -1239,12 +1238,12 @@ C(\boldsymbol{X},\boldsymbol{\theta})=\frac{1}{n}\left\{(\boldsymbol{y}-\boldsym
|
||||
<p>and reordering we have</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\theta}+\frac{n}{2}\lambda sgn(\boldsymbol{\theta})=2\boldsymbol{X}^T\boldsymbol{y}.
|
||||
\boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\theta}+\frac{n}{2}\lambda sgn(\boldsymbol{\theta})=\boldsymbol{X}^T\boldsymbol{y}.
|
||||
\]</div>
|
||||
<p>We can redefine <span class="math notranslate nohighlight">\(\lambda\)</span> to absorb the constant <span class="math notranslate nohighlight">\(n/2\)</span> and we rewrite the last equation as</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\theta}+\lambda sgn(\boldsymbol{\theta})=2\boldsymbol{X}^T\boldsymbol{y}.
|
||||
\boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\theta}+\lambda sgn(\boldsymbol{\theta})=\boldsymbol{X}^T\boldsymbol{y}.
|
||||
\]</div>
|
||||
<p>This equation does not lead to a nice analytical equation as in either Ridge regression or ordinary least squares. This equation can however be solved by using standard convex optimization algorithms.We will discuss how to code the above methods using gradient descent methods.</p>
|
||||
</section>
|
||||
@@ -1643,31 +1642,6 @@ plt.show()
|
||||
</div>
|
||||
</div>
|
||||
</section>
|
||||
<section id="and-a-corresponding-example-using-scikit-learn">
|
||||
<h2>And a corresponding example using <strong>scikit-learn</strong><a class="headerlink" href="#and-a-corresponding-example-using-scikit-learn" title="Link to this heading">#</a></h2>
|
||||
<div class="cell docutils container">
|
||||
<div class="cell_input docutils container">
|
||||
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span># Importing various packages
|
||||
from random import random, seed
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from sklearn.linear_model import SGDRegressor
|
||||
|
||||
n = 100
|
||||
x = 2*np.random.rand(n,1)
|
||||
y = 4+3*x+np.random.randn(n,1)
|
||||
|
||||
X = np.c_[np.ones((n,1)), x]
|
||||
theta_linreg = np.linalg.inv(X.T @ X) @ (X.T @ y)
|
||||
print(theta_linreg)
|
||||
sgdreg = SGDRegressor(max_iter = 50, penalty=None, eta0=0.1)
|
||||
sgdreg.fit(x,y.ravel())
|
||||
print(sgdreg.intercept_, sgdreg.coef_)
|
||||
</pre></div>
|
||||
</div>
|
||||
</div>
|
||||
</div>
|
||||
</section>
|
||||
<section id="gradient-descent-and-ridge">
|
||||
<h2>Gradient descent and Ridge<a class="headerlink" href="#gradient-descent-and-ridge" title="Link to this heading">#</a></h2>
|
||||
<p>We have also discussed Ridge regression where the loss function contains a regularized term given by the <span class="math notranslate nohighlight">\(L_2\)</span> norm of <span class="math notranslate nohighlight">\(\theta\)</span>,</p>
|
||||
@@ -2488,7 +2462,6 @@ plt.show()
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-hessian-matrix">The Hessian matrix</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#simple-program">Simple program</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#id1">Gradient Descent Example</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#and-a-corresponding-example-using-scikit-learn">And a corresponding example using <strong>scikit-learn</strong></a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#gradient-descent-and-ridge">Gradient descent and Ridge</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#the-hessian-matrix-for-ridge-regression">The Hessian matrix for Ridge Regression</a></li>
|
||||
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#program-example-for-gradient-descent-with-ridge-regression">Program example for gradient descent with Ridge Regression</a></li>
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
+365
-403
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,294 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "2c329a04",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
"source": [
|
||||
"<!-- HTML file automatically generated from DocOnce source (https://github.com/doconce/doconce/)\n",
|
||||
"doconce format html exercisesweek37.do.txt -->\n",
|
||||
"<!-- dom:TITLE: Week 37: Linear Regression and Gradient descent -->"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "432a2e5c",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
"source": [
|
||||
"# Week 37: Linear Regression and Gradient descent\n",
|
||||
"**Morten Hjorth-Jensen**, Department of Physics, University of Oslo, Norway\n",
|
||||
"\n",
|
||||
"Date: **September 1-5, 2025**"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "8352513e",
|
||||
"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": "d7eebc41",
|
||||
"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": "1640310e",
|
||||
"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": "60bce352",
|
||||
"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": "75322f31",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
"source": [
|
||||
"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 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": "4b3f945a",
|
||||
"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 = y.mean()\n",
|
||||
"y_centered = y - y_mean"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "12aec8b3",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
"source": [
|
||||
"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": "code",
|
||||
"execution_count": 3,
|
||||
"id": "ffdecd03",
|
||||
"metadata": {
|
||||
"collapsed": false,
|
||||
"editable": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Set regularization parameter\n",
|
||||
"lam = 1.0\n",
|
||||
"\n",
|
||||
"# Closed-form Ridge solution: w = (X^T X + lam * I)^{-1} X^T y\n",
|
||||
"I = np.eye(n_features)\n",
|
||||
"w_closed_form = np.linalg.inv(X_norm.T.dot(X_norm) + lam * I).dot(X_norm.T).dot(y_centered)\n",
|
||||
"\n",
|
||||
"print(\"Closed-form Ridge coefficients:\", w_closed_form)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "33ff39a5",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
"source": [
|
||||
"This computes the ridge 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.\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",
|
||||
"function defined above. The gradient of the ridge cost with respect to\n",
|
||||
"the weight vector $w$ is:\n",
|
||||
"\n",
|
||||
"Below is the code for gradient descent implementation of ridge:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"id": "5b72dc82",
|
||||
"metadata": {
|
||||
"collapsed": false,
|
||||
"editable": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Gradient descent parameters\n",
|
||||
"alpha = 0.1\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 # shape (m,)\n",
|
||||
" # Compute cost (MSE + regularization) for monitoring\n",
|
||||
" cost = (1/(2*m)) * np.dot(error, error) + (lam/(2*m)) * np.dot(theta, theta)\n",
|
||||
" cost_history[t] = cost\n",
|
||||
" # Compute gradient\n",
|
||||
" grad = (1/m) * (X_norm.T.dot(error) + lam * theta)\n",
|
||||
" # Update weights\n",
|
||||
" theta = theta - alpha * grad\n",
|
||||
"\n",
|
||||
"# After the loop, theta contains the fitted coefficients\n",
|
||||
"theta_gd = theta\n",
|
||||
"print(\"Gradient Descent Ridge coefficients:\", theta_gd)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "1fc5b8d0",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
"source": [
|
||||
"Let us confirm that the two approaches (closed-form and gradient\n",
|
||||
"descent) give similar results, and then evaluate the model. First,\n",
|
||||
"compare the learned coefficients to the true coefficients:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"id": "a510e4e6",
|
||||
"metadata": {
|
||||
"collapsed": false,
|
||||
"editable": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"print(\"True coefficients:\", theta_true)\n",
|
||||
"print(\"Closed-form learned coefficients:\", theta_closed_form)\n",
|
||||
"print(\"Gradient descent learned coefficients:\", theta_gd)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "07081487",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
"source": [
|
||||
"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
|
||||
}
|
||||
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,357 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "2c329a04",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
"source": [
|
||||
"<!-- HTML file automatically generated from DocOnce source (https://github.com/doconce/doconce/)\n",
|
||||
"doconce format html exercisesweek37.do.txt -->\n",
|
||||
"<!-- dom:TITLE: Week 37: Linear Regression and Gradient descent -->"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "432a2e5c",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
"source": [
|
||||
"# Week 37: Linear Regression and Gradient descent\n",
|
||||
"**Morten Hjorth-Jensen**, Department of Physics, University of Oslo, Norway\n",
|
||||
"\n",
|
||||
"Date: **September 1-5, 2025**"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "8352513e",
|
||||
"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": 2,
|
||||
"id": "d7eebc41",
|
||||
"metadata": {
|
||||
"collapsed": false,
|
||||
"editable": true,
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"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 @ theta_true + noise"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "1640310e",
|
||||
"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": "60bce352",
|
||||
"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": "75322f31",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
"source": [
|
||||
"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 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": 3,
|
||||
"id": "4b3f945a",
|
||||
"metadata": {
|
||||
"collapsed": false,
|
||||
"editable": true,
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"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 = y.mean()\n",
|
||||
"y_centered = y - y_mean"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "12aec8b3",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
"source": [
|
||||
"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": "code",
|
||||
"execution_count": 7,
|
||||
"id": "ffdecd03",
|
||||
"metadata": {
|
||||
"collapsed": false,
|
||||
"editable": true,
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Closed-form Ridge coefficients: [ 4.98041921 -2.85683951 -0.01549982 0.15658961 -0.07810233 -0.05094777\n",
|
||||
" 1.75377987 0.00679113 0.0343071 -0.04883458]\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# Set regularization parameter\n",
|
||||
"lam = 1.0\n",
|
||||
"\n",
|
||||
"# Closed-form Ridge solution: w = (X^T X + lam * I)^{-1} X^T y\n",
|
||||
"I = np.eye(n_features)\n",
|
||||
"theta_closed_form = np.linalg.inv(X_norm.T.dot(X_norm) + lam * I).dot(X_norm.T).dot(y_centered)\n",
|
||||
"\n",
|
||||
"print(\"Closed-form Ridge coefficients:\", theta_closed_form)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "33ff39a5",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
"source": [
|
||||
"This computes the ridge 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.\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",
|
||||
"function defined above. The gradient of the ridge cost with respect to\n",
|
||||
"the weight vector $w$ is:\n",
|
||||
"\n",
|
||||
"Below is the code for gradient descent implementation of ridge:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"id": "5b72dc82",
|
||||
"metadata": {
|
||||
"collapsed": false,
|
||||
"editable": true,
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Gradient Descent Ridge coefficients: [ 4.98041921 -2.85683951 -0.01549982 0.15658961 -0.07810233 -0.05094777\n",
|
||||
" 1.75377987 0.00679113 0.0343071 -0.04883458]\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# Gradient descent parameters\n",
|
||||
"alpha = 0.1\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 # shape (m,)\n",
|
||||
" # Compute cost (MSE + regularization) for monitoring\n",
|
||||
" cost = (1/(2*m)) * np.dot(error, error) + (lam/(2*m)) * np.dot(theta, theta)\n",
|
||||
" cost_history[t] = cost\n",
|
||||
" # Compute gradient\n",
|
||||
" grad = (1/m) * (X_norm.T.dot(error) + lam * theta)\n",
|
||||
" # Update weights\n",
|
||||
" theta = theta - alpha * grad\n",
|
||||
"\n",
|
||||
"# After the loop, theta contains the fitted coefficients\n",
|
||||
"theta_gd = theta\n",
|
||||
"print(\"Gradient Descent Ridge coefficients:\", theta_gd)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "1fc5b8d0",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
"source": [
|
||||
"Let us confirm that the two approaches (closed-form and gradient\n",
|
||||
"descent) give similar results, and then evaluate the model. First,\n",
|
||||
"compare the learned coefficients to the true coefficients:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"id": "a510e4e6",
|
||||
"metadata": {
|
||||
"collapsed": false,
|
||||
"editable": true,
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"True coefficients: [ 5. -3. 0. 0. 0. 0. 2. 0. 0. 0.]\n",
|
||||
"Closed-form learned coefficients: [ 4.98041921 -2.85683951 -0.01549982 0.15658961 -0.07810233 -0.05094777\n",
|
||||
" 1.75377987 0.00679113 0.0343071 -0.04883458]\n",
|
||||
"Gradient descent learned coefficients: [ 4.98041921 -2.85683951 -0.01549982 0.15658961 -0.07810233 -0.05094777\n",
|
||||
" 1.75377987 0.00679113 0.0343071 -0.04883458]\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"print(\"True coefficients:\", theta_true)\n",
|
||||
"print(\"Closed-form learned coefficients:\", theta_closed_form)\n",
|
||||
"print(\"Gradient descent learned coefficients:\", theta_gd)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "07081487",
|
||||
"metadata": {
|
||||
"editable": true
|
||||
},
|
||||
"source": [
|
||||
"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": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 3 (ipykernel)",
|
||||
"language": "python",
|
||||
"name": "python3"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 3
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.9.15"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 5
|
||||
}
|
||||
Reference in New Issue
Block a user