Multiple Linear Regression

1. From a Line to a Hyperplane

In simple linear regression, our model was a line described by the equation:
$\hat{y}={\beta }_0+{\beta }_1{x}_1$
Here, ${\beta }_0$ is the intercept and ${\beta }_1$ is the coefficient (slope) for our single feature, $x_1$.

To incorporate more features, we simply extend this equation by adding more terms, one for each new feature:
$\hat{y}={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+{\beta }_3{x}_3+\cdots +{\beta }_p{x}_p$

• $\hat{y}$ is still our predicted value.
• $x_1,x_2,\dots ,x_p$ are the $p$ different features.
• ${\beta }_0$ is the intercept.
• ${\beta }_1,{\beta }_2,\dots ,{\beta }_p$ are the coefficients for each feature. Each ${\beta }_j$ represents the change in $\hat{y}$ for a one-unit change in $x_j$, assuming all other features are held constant.

Visually, we're no longer fitting a 2D line. With two features ($x_1,x_2$), we're fitting a 3D plane. With more than two features, we're fitting a hyperplane, which is just the general term for a flat surface in a multi-dimensional space.

To get a clearer picture of this transition, please start by reading the first two short sections of the following guide.

2. The Power of Matrix Notation

As you can see, the equation gets long and clumsy with many features. Your computer science background will appreciate that we can represent this much more elegantly using linear algebra.

The entire system of equations for all our data points can be expressed in a single, compact form:
$\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$

Let's break down each component:

• $\mathbf{y}$ is a vector containing the actual target values for all $n$ observations in our dataset.
• $\mathbf{\beta }$ is a vector containing all our model's coefficients (${\beta }_0,{\beta }_1,\dots ,{\beta }_p$).
• $\mathbf{ϵ}$ is a vector containing the error terms (residuals) for each observation.
• $\mathbf{X}$ is the design matrix. This is a crucial concept. It's an $n\times (p+1)$ matrix where each row is an observation and each column is a feature. We add a special first column of all ones. This "dummy" column corresponds to the intercept term ${\beta }_0$, ensuring it gets estimated correctly.

The video below provides an excellent visual breakdown of how the individual data points are organized into these vectors and the design matrix.

3. Finding the Optimal Coefficients

Just like in simple linear regression, our goal is to find the coefficient vector $\mathbf{\beta }$ that minimizes the sum of squared errors (SSE).
$SSE={\sum }_{i=1}^n(y_i-{\hat{y}}_i{)}^2$
In matrix form, this is written as:
$SSE=(\mathbf{y}-\mathbf{X}\mathbf{\beta }{)}^T(\mathbf{y}-\mathbf{X}\mathbf{\beta })$

There are two primary methods to find the $\mathbf{\beta }$ that minimizes this SSE.

Method 1: The Normal Equation (Analytical Solution)

This is the multi-feature equivalent of the OLS formulas we derived in the last lesson. By using matrix calculus to take the derivative of the SSE with respect to $\mathbf{\beta }$ and setting it to zero, we arrive at a beautiful closed-form solution called the Normal Equation:

$\hat{\mathbf{\beta }}=({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T\mathbf{y}$

• $\hat{\mathbf{\beta }}$ is our estimated coefficient vector.
• ${\mathbf{X}}^T$ is the transpose of the design matrix $\mathbf{X}$.
• $({\mathbf{X}}^T\mathbf{X}{)}^{-1}$ is the inverse of the matrix product ${\mathbf{X}}^T\mathbf{X}$.

This equation gives us the exact optimal coefficients in one direct calculation!

To see how this equation is derived and used, continue with the video from Boer Commander.

The main advantage of the Normal Equation is its directness. However, calculating the matrix inverse $({\mathbf{X}}^T\mathbf{X}{)}^{-1}$ is computationally very expensive, especially for datasets with a large number of features. Its complexity is roughly $O(p^3)$, where $p$ is the number of features. This makes it impractical for problems with tens of thousands of features.

Method 2: Gradient Descent (Iterative Solution)

When the Normal Equation is too slow, we turn to gradient descent. As you'll recall from our Calculus module, this is an iterative optimization algorithm that finds the minimum of a function.

The process is:

1. Start with some initial values for the coefficients in $\mathbf{\beta }$ (e.g., all zeros).
2. Calculate the cost (SSE) for the current coefficients.
3. Calculate the gradient (the partial derivative of the cost function with respect to each coefficient). This tells us the direction of steepest ascent.
4. Update each coefficient by taking a small step in the opposite direction of the gradient. The size of this step is controlled by the learning rate, $\alpha$.
5. Repeat steps 2-4 until the cost function converges to a minimum.

The update rule for each coefficient ${\beta }_j$ in the vector $\mathbf{\beta }$ is:
${\beta }_j:={\beta }_j-\alpha \frac{\partial }{\partial {\beta }_j}J(\mathbf{\beta })$
where $J(\mathbf{\beta })$ is our cost function.

The following video provides an excellent and in-depth explanation of linear regression with multiple features and how gradient descent is used to optimize the parameters.

In practice, for very large datasets, a variant called Stochastic Gradient Descent (SGD) is often used. Instead of calculating the gradient using the entire dataset in each step (which is slow), SGD updates the parameters using just one training example at a time. This makes progress much faster, though the path to the minimum is "noisier." The Stanford video mentions this toward the end (starting at 44:19), which you can explore if you're curious, but it's not essential for this lesson.

4. Implementation in Python

Now, let's put this theory into code. We'll use the excellent guide from mbrenndoerfer.com to see how to implement MLR.

First, we'll see how to do it with scikit-learn, which is the standard practical approach. Then, for a deeper understanding, we will implement the Normal Equation from scratch using NumPy.

Here's the key code block from that article for implementing the Normal Equation with NumPy. It's a direct translation of the math into code.

import numpy as np

def multiple_linear_regression_normal_equation(X, y):
    """
    Fits multiple linear regression using the closed-form OLS solution (Normal Equation).
    
    Parameters:
    X: Feature matrix (n_samples, n_features)
    y: Target vector (n_samples,)
    
    Returns:
    coefficients: Array of coefficients [intercept, feature1, feature2, ...]
    """
    # Add a column of ones to X for the intercept term
    X_with_intercept = np.c_[np.ones(X.shape[0]), X]
    
    # Compute the OLS solution: β = (X'X)^(-1) X'y
    XtX = X_with_intercept.T @ X_with_intercept
    Xty = X_with_intercept.T @ y
    
    # Use np.linalg.solve for better numerical stability and efficiency
    # This solves the equation XtX * β = Xty for β
    coefficients = np.linalg.solve(XtX, Xty)
    
    return coefficients

# --- Let's test it on a small dataset ---
# Features: [size in sq ft, age in years]
X_np = np.array([[2000, 1], [1500, 5], [2500, 2]])
# Target: price in thousands
y_np = np.array([400, 300, 500])

# Get the coefficients
coeffs = multiple_linear_regression_normal_equation(X_np, y_np)
intercept, slope_size, slope_age = coeffs[0], coeffs[1], coeffs[2]

print(f"Calculated coefficients using Normal Equation:")
print(f"Intercept (β0): {intercept:.4f}")
print(f"Size coefficient (β1): {slope_size:.4f}")
print(f"Age coefficient (β2): {slope_age:.4f}")

# Compare with scikit-learn
from sklearn.linear_model import LinearRegression

model = LinearRegression()
model.fit(X_np, y_np)

print("\nScikit-learn's coefficients:")
print(f"Intercept (β0): {model.intercept_:.4f}")
print(f"Coefficients (β1, β2): {model.coef_}")

Running this code confirms that our from-scratch implementation of the Normal Equation yields the same results as scikit-learn, proving our understanding of the underlying mechanics.

Conclusion

In this lesson, you've made a significant leap from simple to multiple linear regression, unlocking the ability to model more complex, real-world relationships.

Key Takeaways:

• Multiple Linear Regression extends the simple model by adding a coefficient for each additional feature, fitting a hyperplane to the data.
• The equation can be expressed compactly using matrix notation: $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$. The design matrix $\mathbf{X}$, with its leading column of ones, is a key component.
• There are two main ways to find the optimal coefficients $\mathbf{\beta }$ that minimize SSE:
  • The Normal Equation: A direct, analytical solution that is elegant but computationally expensive for many features.
  • Gradient Descent: An iterative optimization approach that is more scalable and widely used in modern machine learning.

Preview of the Next Lesson:
You might think linear regression can only model linear relationships, but that's not quite true! In our next lesson, we will explore Polynomial Regression. You'll learn a clever technique to use the linear regression framework we've just studied to fit non-linear, curved data by creating new features from our existing ones.