> ## Documentation Index
> Fetch the complete documentation index at: https://resources.devweekends.com/llms.txt
> Use this file to discover all available pages before exploring further.
# Linear Regression
> Your first complete ML algorithm - fitting lines through points
# Linear Regression
## From Intuition to Implementation
In the previous modules, you learned:
1. ML is about finding patterns in data
2. We measure "wrongness" with a loss function
3. Gradient descent minimizes the loss
Now let's put it all together into a complete, professional algorithm.
***
## The Real-World Setup
**Your boss asks**: "Can you predict how much revenue we'll make based on our advertising spend?"
You have historical data:
```python theme={null}
import numpy as np
import matplotlib.pyplot as plt
# Historical data: [TV ads ($k), Radio ads ($k), Newspaper ads ($k)] -> Revenue ($k)
advertising = np.array([
[230, 37, 69],
[44, 39, 45],
[17, 45, 69],
[151, 41, 58],
[180, 10, 58],
[8, 49, 75],
[57, 32, 23],
[120, 19, 11],
[8, 35, 65],
[199, 3, 70],
])
revenue = np.array([22.1, 10.4, 9.3, 18.5, 12.9, 7.2, 11.8, 13.2, 4.8, 10.6])
```
**Question**: If we spend $100k on TV, $40k on Radio, and \$30k on Newspaper, what revenue can we expect?
***
## The Linear Regression Model
We assume revenue is a weighted combination of ad spends:
$$
\text{revenue} = w_0 + w_1 \cdot \text{TV} + w_2 \cdot \text{Radio} + w_3 \cdot \text{Newspaper}
$$
Or in matrix notation (see [Matrix Operations](/courses/math-for-ml-linear-algebra/03-matrices)):
$$
\hat{y} = X \cdot w
$$
Where:
* $X$ is our feature matrix (with a column of 1s for the bias term)
* $w$ is our weight vector
* $\hat{y}$ is our predictions
***
## Step-by-Step Implementation
### Step 1: Prepare the Data
```python theme={null}
# Add a column of 1s for the bias (intercept) term
X = np.column_stack([np.ones(len(advertising)), advertising])
y = revenue
print("X shape:", X.shape) # (10, 4) - 10 samples, 4 features (including bias)
print("y shape:", y.shape) # (10,) - 10 target values
```
### Step 2: Define the Model
```python theme={null}
def predict(X, w):
"""Linear model: y_hat = X @ w"""
return X @ w
```
### Step 3: Define the Loss Function
```python theme={null}
def mean_squared_error(y_true, y_pred):
"""Average of squared differences."""
return np.mean((y_true - y_pred) ** 2)
def compute_loss(X, y, w):
"""Compute MSE for given weights."""
predictions = predict(X, w)
return mean_squared_error(y, predictions)
```
### Step 4: Compute the Gradient
```python theme={null}
def compute_gradient(X, y, w):
"""
Gradient of MSE with respect to weights.
d(MSE)/dw = (2/n) * X.T @ (X @ w - y)
"""
n = len(y)
predictions = predict(X, w)
errors = predictions - y
gradient = (2 / n) * X.T @ errors
return gradient
```
### Step 5: Gradient Descent
```python theme={null}
def train_linear_regression(X, y, learning_rate=0.0001, num_epochs=1000):
"""
Train a linear regression model using gradient descent.
Args:
X: Feature matrix (n_samples, n_features)
y: Target values (n_samples,)
learning_rate: Step size for gradient descent
num_epochs: Number of training iterations
Returns:
Trained weights
"""
# Initialize weights to zero
w = np.zeros(X.shape[1])
# Track loss history for plotting
losses = []
for epoch in range(num_epochs):
# Compute current loss
loss = compute_loss(X, y, w)
losses.append(loss)
# Compute gradient
gradient = compute_gradient(X, y, w)
# Update weights
w = w - learning_rate * gradient
if epoch % 100 == 0:
print(f"Epoch {epoch}: Loss = {loss:.4f}")
return w, losses
```
### Step 6: Train the Model
```python theme={null}
# Normalize features for better convergence
X_normalized = X.copy()
X_normalized[:, 1:] = (X[:, 1:] - X[:, 1:].mean(axis=0)) / X[:, 1:].std(axis=0)
# Train
weights, loss_history = train_linear_regression(X_normalized, y, learning_rate=0.1, num_epochs=500)
print("\nFinal weights:")
print(f" Bias: {weights[0]:.4f}")
print(f" TV: {weights[1]:.4f}")
print(f" Radio: {weights[2]:.4f}")
print(f" Newspaper: {weights[3]:.4f}")
```
***
## Using scikit-learn (The Professional Way)
In practice, we use libraries that handle all the details -- the normal equation, numerical stability, edge cases. Understanding the math (above) is for your brain; scikit-learn is for your production code.
```python theme={null}
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, r2_score
# Split data into training and test sets.
# Why 80/20? It's a common default. With very little data, consider
# cross-validation instead (Module 7) to use every row for both
# training and evaluation.
X_train, X_test, y_train, y_test = train_test_split(
advertising, revenue, test_size=0.2, random_state=42
)
# Create and train the model.
# Under the hood, sklearn uses the normal equation (not gradient descent)
# for LinearRegression -- it's faster for small datasets.
model = LinearRegression()
model.fit(X_train, y_train)
# Make predictions on data the model has NEVER seen
y_pred = model.predict(X_test)
# Evaluate: R-squared tells you what fraction of the variance
# in revenue is explained by ad spend. 0.8 means "80% explained."
mse = mean_squared_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)
print(f"Mean Squared Error: {mse:.4f}")
print(f"R-squared Score: {r2:.4f}")
# The coefficients are the "weights" we've been learning about.
# Each one tells you: "holding everything else constant,
# how much does a $1K increase in this ad channel change revenue?"
print("\nLearned Coefficients:")
print(f" Intercept: {model.intercept_:.4f}")
print(f" TV: {model.coef_[0]:.4f}")
print(f" Radio: {model.coef_[1]:.4f}")
print(f" Newspaper: {model.coef_[2]:.4f}")
```
***
## Interpreting the Results
The coefficients tell a story:
```
TV coefficient: 0.046
Radio coefficient: 0.188
Newspaper coefficient: 0.003
```
**Insights**:
* Every $1k spent on **Radio** returns ~$188 in revenue (best ROI!)
* Every $1k on **TV** returns ~$46
* **Newspaper** ads have almost no impact
**Business decision**: Shift budget from newspaper to radio!
***
## The Closed-Form Solution
For linear regression, there's actually a formula that gives the optimal weights directly, without gradient descent:
$$
w = (X^T X)^{-1} X^T y
$$
This is called the **Normal Equation**. It comes from calculus - setting the gradient to zero and solving. Think of it as the "just give me the answer" approach versus gradient descent's "let me walk there step by step."
```python theme={null}
def linear_regression_closed_form(X, y):
"""
Compute optimal weights using the normal equation.
This solves the system of equations directly -- no iterations,
no learning rate to tune. The trade-off? It requires computing
a matrix inverse, which is expensive for many features.
"""
# w = (X^T X)^(-1) X^T y
XtX = X.T @ X # Shape: (features, features) -- square matrix
XtX_inv = np.linalg.inv(XtX) # O(n^3) -- the expensive step
Xty = X.T @ y # Shape: (features,) -- project targets onto features
w = XtX_inv @ Xty
return w
# This gives the same answer as gradient descent!
w_closed = linear_regression_closed_form(X, y)
print("Closed-form solution:", w_closed)
```
**When to use which?**
| Approach | Best When | Why |
| -------------------- | ---------------------------------------- | -------------------------------------------- |
| **Normal Equation** | \< 10,000 features | Exact solution, no hyperparameters to tune |
| **Gradient Descent** | > 10,000 features or very large datasets | O(n) per step vs O(n^3) for matrix inversion |
In practice, scikit-learn's `LinearRegression` automatically picks the best solver for your data size. You rarely need to worry about this choice -- but understanding it helps you debug slow training times.
***
## Real-World Example: House Price Prediction
Let's build a proper house price predictor using real data:
```python theme={null}
from sklearn.datasets import fetch_california_housing
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_squared_error, r2_score
# Load California housing data
housing = fetch_california_housing()
X, y = housing.data, housing.target
print("Features:", housing.feature_names)
print("X shape:", X.shape) # (20640, 8) - 20k houses, 8 features
# Split data
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.2, random_state=42
)
# Scale features (important for many algorithms!)
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
# Train model
model = LinearRegression()
model.fit(X_train_scaled, y_train)
# Evaluate
y_pred = model.predict(X_test_scaled)
print(f"\nR-squared: {r2_score(y_test, y_pred):.4f}")
print(f"RMSE: ${np.sqrt(mean_squared_error(y_test, y_pred)) * 100000:.2f}")
# Feature importance
print("\nFeature Importance (coefficients):")
for name, coef in zip(housing.feature_names, model.coef_):
print(f" {name}: {coef:.4f}")
```
***
## Common Pitfalls and Solutions
### Pitfall 1: Features on Different Scales
```python theme={null}
# BAD: Sqft ranges 500-5000, bedrooms 1-6
# The model will be biased toward sqft because the
# coefficient for sqft will be tiny (e.g., 0.15) while
# the bedroom coefficient will be huge (e.g., 25000).
# This doesn't affect predictions, but it makes coefficients
# impossible to compare for feature importance.
# GOOD: Standardize features so each has mean=0, std=1.
# Now a 1-unit change in any feature means "1 standard deviation,"
# making coefficients directly comparable.
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
```
### Pitfall 2: Multicollinearity
When features are highly correlated, coefficients become unstable. Imagine trying to figure out whether it's the coffee or the sugar making your drink sweet -- when they always appear together, it's hard to tell who deserves the credit.
```python theme={null}
# Check correlations -- look for values above 0.8 or below -0.8
import pandas as pd
df = pd.DataFrame(X, columns=feature_names)
print(df.corr())
# Solution: Use Ridge regression (L2 regularization) to stabilize
# coefficients, or drop one of the correlated features.
# Ridge adds a penalty that discourages any single weight from
# getting too large, which naturally handles collinearity.
from sklearn.linear_model import Ridge
model = Ridge(alpha=1.0) # alpha controls regularization strength
```
### Pitfall 3: Overfitting
When the model memorizes training data but fails on new data. This is like studying only the practice exam and then bombing the real test because the questions are slightly different.
```python theme={null}
# Always evaluate on held-out test data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
model.fit(X_train, y_train)
print("Train R2:", model.score(X_train, y_train))
print("Test R2:", model.score(X_test, y_test))
# If train >> test, you're overfitting!
# A 5-10% gap is normal. A 30%+ gap means trouble.
```
**Practical rule of thumb**: Linear regression rarely overfits unless you have far more features than samples. If you have 50 features and 100 rows, consider Ridge or Lasso regression (Module 13) to keep things under control.
***
## The Complete Linear Regression Workflow
```python theme={null}
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.metrics import mean_squared_error, r2_score
def linear_regression_pipeline(X, y, feature_names=None):
"""
Complete linear regression workflow.
"""
# 1. Split data
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.2, random_state=42
)
# 2. Scale features
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
# 3. Train model
model = LinearRegression()
model.fit(X_train_scaled, y_train)
# 4. Evaluate
y_pred_train = model.predict(X_train_scaled)
y_pred_test = model.predict(X_test_scaled)
print("=== Model Evaluation ===")
print(f"Train R2: {r2_score(y_train, y_pred_train):.4f}")
print(f"Test R2: {r2_score(y_test, y_pred_test):.4f}")
print(f"Test RMSE: {np.sqrt(mean_squared_error(y_test, y_pred_test)):.4f}")
# 5. Cross-validation
cv_scores = cross_val_score(model, X_train_scaled, y_train, cv=5)
print(f"\nCross-validation R2: {cv_scores.mean():.4f} (+/- {cv_scores.std():.4f})")
# 6. Feature importance
if feature_names:
print("\n=== Feature Importance ===")
importance = pd.DataFrame({
'feature': feature_names,
'coefficient': model.coef_
}).sort_values('coefficient', key=abs, ascending=False)
print(importance.to_string(index=False))
return model, scaler
```
***
## Key Takeaways
y = w0 + w1*x1 + w2*x2 + ...
Measures average squared error
Normalize for better training
Always hold out some data
***
## 🚀 Mini Projects
Build a salary prediction model
Real estate price predictor with feature engineering
Model comparison and selection pipeline
**Project 1: Salary Prediction Model** - Complete regression workflow
**Objective**: Predict software engineer salaries based on experience, skills, and location.
```python theme={null}
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.preprocessing import StandardScaler
# Generate realistic salary data
np.random.seed(42)
n = 200
years_exp = np.random.uniform(0, 15, n)
num_skills = np.random.randint(2, 12, n)
is_sf = np.random.choice([0, 1], n, p=[0.7, 0.3])
has_masters = np.random.choice([0, 1], n, p=[0.75, 0.25])
# True salary model with noise
salary = (50000 +
7000 * years_exp +
3000 * num_skills +
25000 * is_sf +
15000 * has_masters +
np.random.normal(0, 10000, n))
# Create feature matrix
X = np.column_stack([years_exp, num_skills, is_sf, has_masters])
y = salary
# Split data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# Scale features
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
# Train model
model = LinearRegression()
model.fit(X_train_scaled, y_train)
# Evaluate
y_pred = model.predict(X_test_scaled)
rmse = np.sqrt(mean_squared_error(y_test, y_pred))
r2 = r2_score(y_test, y_pred)
print("=== Model Performance ===")
print(f"R² Score: {r2:.4f}")
print(f"RMSE: ${rmse:,.0f}")
# Feature importance
feature_names = ['Years Exp', 'Num Skills', 'SF Location', 'Masters Degree']
print("\n=== Feature Importance ===")
for name, coef in zip(feature_names, model.coef_):
print(f"{name}: {coef:+,.0f}")
# Predictions for specific profiles
print("\n=== Salary Predictions ===")
profiles = [
[2, 5, 0, 0], # Junior in non-SF
[5, 7, 1, 0], # Mid-level in SF
[10, 10, 1, 1], # Senior in SF with Masters
]
for profile in profiles:
scaled = scaler.transform([profile])
pred = model.predict(scaled)[0]
print(f"Exp={profile[0]}yr, Skills={profile[1]}, SF={profile[2]}, Masters={profile[3]}")
print(f" → Predicted: ${pred:,.0f}\n")
```
**Project 2: Real Estate Price Predictor** - Feature engineering
**Objective**: Build a house price predictor with feature engineering.
```python theme={null}
import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.preprocessing import StandardScaler, PolynomialFeatures
from sklearn.metrics import r2_score
# Generate housing data
np.random.seed(42)
n = 300
sqft = np.random.normal(2000, 500, n)
bedrooms = np.random.choice([2, 3, 4, 5], n, p=[0.1, 0.4, 0.4, 0.1])
bathrooms = np.random.choice([1, 2, 3], n, p=[0.2, 0.5, 0.3])
age = np.random.exponential(20, n)
lot_size = np.random.normal(8000, 2000, n)
# Non-linear price relationship
price = (100000 +
150 * sqft +
20000 * bedrooms +
15000 * bathrooms -
1500 * age -
0.001 * age**2 * sqft + # Interaction: old big houses depreciate more
5 * lot_size +
np.random.normal(0, 30000, n))
# Create DataFrame
df = pd.DataFrame({
'sqft': sqft, 'bedrooms': bedrooms, 'bathrooms': bathrooms,
'age': age, 'lot_size': lot_size, 'price': price
})
# Feature Engineering
df['price_per_sqft'] = df['price'] / df['sqft'] # Target leak - don't use!
df['sqft_per_bedroom'] = df['sqft'] / df['bedrooms']
df['bath_to_bed_ratio'] = df['bathrooms'] / df['bedrooms']
df['is_new'] = (df['age'] < 5).astype(int)
df['is_large'] = (df['sqft'] > 2500).astype(int)
# Features to use
feature_cols = ['sqft', 'bedrooms', 'bathrooms', 'age', 'lot_size',
'sqft_per_bedroom', 'bath_to_bed_ratio', 'is_new', 'is_large']
X = df[feature_cols].values
y = df['price'].values
# Compare models
print("=== Model Comparison ===\n")
# Model 1: Basic features only
X_basic = df[['sqft', 'bedrooms', 'bathrooms', 'age']].values
X_train, X_test, y_train, y_test = train_test_split(X_basic, y, test_size=0.2, random_state=42)
model_basic = LinearRegression().fit(X_train, y_train)
r2_basic = r2_score(y_test, model_basic.predict(X_test))
print(f"Basic features: R² = {r2_basic:.4f}")
# Model 2: With engineered features
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
model_eng = LinearRegression().fit(X_train, y_train)
r2_eng = r2_score(y_test, model_eng.predict(X_test))
print(f"Engineered features: R² = {r2_eng:.4f}")
# Model 3: Polynomial features
poly = PolynomialFeatures(degree=2, include_bias=False)
X_poly = poly.fit_transform(df[['sqft', 'age']].values)
X_train, X_test, y_train, y_test = train_test_split(X_poly, y, test_size=0.2, random_state=42)
model_poly = LinearRegression().fit(X_train, y_train)
r2_poly = r2_score(y_test, model_poly.predict(X_test))
print(f"Polynomial features: R² = {r2_poly:.4f}")
print(f"\nImprovement from engineering: {(r2_eng - r2_basic) * 100:.1f}%")
```
**Project 3: Model Comparison Pipeline** - End-to-end ML workflow
**Objective**: Compare multiple models and select the best one.
```python theme={null}
import numpy as np
from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet
from sklearn.preprocessing import StandardScaler, PolynomialFeatures
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.pipeline import Pipeline
# Generate data with some noise
np.random.seed(42)
n = 500
X = np.random.randn(n, 5)
# True relationship with only 3 relevant features
y = 3*X[:, 0] + 2*X[:, 1] - 1.5*X[:, 2] + np.random.randn(n) * 0.5
# Split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# Define models to compare
models = {
'Linear Regression': LinearRegression(),
'Ridge (α=1.0)': Ridge(alpha=1.0),
'Ridge (α=10.0)': Ridge(alpha=10.0),
'Lasso (α=0.1)': Lasso(alpha=0.1),
'Lasso (α=1.0)': Lasso(alpha=1.0),
'ElasticNet': ElasticNet(alpha=0.1, l1_ratio=0.5),
}
print("=== Model Comparison ===\n")
print(f"{'Model':<25} {'Train R²':>10} {'Test R²':>10} {'CV R² (±std)':>18}")
print("-" * 65)
results = []
for name, model in models.items():
# Fit
model.fit(X_train, y_train)
# Scores
train_r2 = r2_score(y_train, model.predict(X_train))
test_r2 = r2_score(y_test, model.predict(X_test))
cv_scores = cross_val_score(model, X_train, y_train, cv=5)
print(f"{name:<25} {train_r2:>10.4f} {test_r2:>10.4f} {cv_scores.mean():>8.4f} (±{cv_scores.std():.4f})")
results.append((name, test_r2, model))
# Best model
best_name, best_r2, best_model = max(results, key=lambda x: x[1])
print(f"\n=== Best Model: {best_name} ===")
print(f"Test R²: {best_r2:.4f}")
# Feature importance from best model
print("\n=== Coefficients ===")
for i, coef in enumerate(best_model.coef_):
marker = "***" if abs(coef) > 0.5 else ""
print(f"Feature {i+1}: {coef:+.4f} {marker}")
# Lasso feature selection
lasso = Lasso(alpha=0.5)
lasso.fit(X_train, y_train)
print("\n=== Lasso Feature Selection ===")
print("Non-zero coefficients indicate important features:")
for i, coef in enumerate(lasso.coef_):
if abs(coef) > 0.01:
print(f" Feature {i+1}: {coef:+.4f}")
```
***
## What's Next?
Linear regression predicts continuous numbers. But what if you want to predict **categories**?
* Is this email spam or not spam?
* Will this customer churn or stay?
* Is this tumor malignant or benign?
That's **classification** - the subject of our next module!
Learn to predict categories with logistic regression and beyond
***
## 🔗 Math → ML Connection Summary
**Where the math you learned powers linear regression:**
| Math Concept | Where It Appears | What It Does |
| ------------------------------------------------------------------------------------- | ----------------------------- | ----------------------------------- |
| **Vectors** ([Module](/courses/math-for-ml-linear-algebra/02-vectors)) | Feature vector $\mathbf{x}$ | Represents each data point |
| **Dot Product** ([Module](/courses/math-for-ml-linear-algebra/02-vectors)) | $\mathbf{w} \cdot \mathbf{x}$ | Core of the linear model prediction |
| **Matrix Multiplication** ([Module](/courses/math-for-ml-linear-algebra/03-matrices)) | $X\mathbf{w} = \hat{y}$ | Batch predictions for all samples |
| **Gradient** ([Module](/courses/math-for-ml-calculus/02-gradients)) | $\nabla_w L$ | Direction to update weights |
| **Chain Rule** ([Module](/courses/math-for-ml-calculus/03-chain-rule)) | Backpropagation | Compute gradients efficiently |
| **Mean/Variance** ([Module](/courses/statistics-for-ml/02-describing-data)) | Feature scaling | Standardization for better training |
| **Normal Distribution** ([Module](/courses/statistics-for-ml/04-distributions)) | Error distribution | Justifies using MSE loss |
**Bottom line**: Linear regression is the intersection of all three math courses. Master this, and neural networks become "just deeper linear regression with nonlinearities."
**Want to understand the theory?** Here's what's happening under the hood:
### Why Gradient Descent Works
The MSE loss is **convex** (bowl-shaped), meaning:
* There's exactly one minimum
* Gradient descent is guaranteed to find it
* Step size (learning rate) affects speed but not destination
$L(w) = \frac{1}{n}\sum_{i=1}^{n}(y_i - \mathbf{w}^T\mathbf{x}_i)^2$
Taking the gradient: $\nabla_w L = -\frac{2}{n}X^T(y - X\mathbf{w})$
Setting to zero gives the normal equation: $\mathbf{w}^* = (X^TX)^{-1}X^Ty$
### The Statistical Interpretation
Under the assumption that errors are normally distributed:
$y = \mathbf{w}^T\mathbf{x} + \epsilon, \quad \epsilon \sim \mathcal{N}(0, \sigma^2)$
**Maximum Likelihood Estimation** (MLE) of $\mathbf{w}$ is equivalent to minimizing MSE!
### Regularization from a Bayesian View
* **Ridge Regression** (L2): Assumes weights have Gaussian prior $\mathbf{w} \sim \mathcal{N}(0, \tau^2I)$
* **Lasso** (L1): Assumes weights have Laplacian prior → promotes sparsity
This connection between regularization and Bayesian priors is why regularization prevents overfitting!
### Recommended Deep-Dive Resources
* *The Elements of Statistical Learning* - Chapter 3 (free online)
* *Pattern Recognition and Machine Learning* by Bishop - Chapter 3
* [Stanford CS229 Linear Regression Notes](https://cs229.stanford.edu/notes2021fall/cs229-notes1.pdf)