Linear Regression and the Normal Equation

Linear Regression models a continuous target $y$ as a linear combination of input features $x$:

$\hat{y} = \theta^T \mathbf{x}$ where $\theta$ is the model parameter vector (including the intercept or bias $\theta_0$), and $\mathbf{x}$ is the observation vector containing the features (with $x_0 = 1$).

Simple vs. Multiple Linear Regression

• Simple Linear Regression: Models the relationship between a single feature $x_1$ and the target $y$. The equation simplifies to finding the best line with a specific slope ($\theta_1$) and intercept ($\theta_0$): $\hat{y} = \theta_0 + \theta_1 x_1$
• Multiple Linear Regression: Extends the concept to multiple features, where each feature gets its own weight (slope) in a multidimensional space.

Polynomial Regression

If the data is more complex than a simple straight line, we can still use linear models to fit nonlinear data. By adding powers of each feature as new features, we can train a linear model on this extended set of features. This is called Polynomial Regression. Despite fitting a curve to the data, it is still considered a linear model because the prediction is still a linear combination of the (now polynomial) features.

Parameter Estimation: Ordinary Least Squares

To train the model, we find parameters $\theta$ that minimize the Mean Squared Error (MSE) over the dataset:

$\text{MSE}(X, h_{\theta}) = \frac{1}{m} \sum_{i=1}^{m} (\theta^T \mathbf{x}^{(i)} - y^{(i)})^2$

The Closed-Form Normal Equation

To minimize the cost function, we solve analytically using the Normal Equation:

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

• $X$: The design matrix of shape $(m, n+1)$ containing all features.
• $y$: The target vector of shape $(m, 1)$.

Multicollinearity and SVD

If the matrix $X^T X$ is singular (non-invertible) due to redundant, highly correlated features, standard inversion fails. Solvers compute the pseudo-inverse $X^+$ using Singular Value Decomposition (SVD) for stability.

Computational Complexity

Computing $(X^T X)^{-1}$ requires inverting an $(n+1) \times (n+1)$ matrix. The computational complexity is between $O(n^{2.4})$ and $O(n^3)$, making the Normal Equation expensive when the number of features is large.

Example: Computing the Normal Equation

The following example demonstrates computing model coefficients using the Normal Equation:

import numpy as np

# Generate synthetic linear data with random Gaussian noise
np.random.seed(42)
X = 2 * np.random.rand(5, 1)
y = 4 + 3 * X + np.random.randn(5, 1)

# Add bias column (x0 = 1) to build design matrix X_b
X_b = np.c_[np.ones((5, 1)), X]

# Solve analytically via normal equation: theta = (X^T * X)^(-1) * X^T * y
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)
print("Estimated coefficients [bias, weight]:")
print(theta_best.flatten())

Estimated coefficients [bias, weight]:
[3.9877443  3.40798606]

Exercise

Test your understanding of the analytical parameter estimation limits:

References & Further Reading

• Géron, A. (2022). Hands-On Machine Learning with Scikit-Learn, Keras, and TensorFlow (3rd ed.). O’Reilly Media.
• Smola, A., & Vishwanathan, S. V. N. (2008). An Introduction to Machine Learning. Cambridge University Press.