Linear models from least squares to quantiles

In the last few posts we built up the machinery of supervised learning from first principles, from nearest neighbors through decision trees and ensemble methods. Now we turn to the oldest and most enduring family in the toolbox: linear models. By the end of this post we will have fitted the major linear regression families to the California Housing dataset, watched our test error fall from a naive baseline of 1.1449 to a polynomial ridge model at 0.6877, and learned why the choice of loss function matters as much as the choice of model.

The through-line for this post is a simple one: every linear model is the same machine, and the only thing that changes is what we tell it to care about. Least squares cares about every error equally. Huber loss cares less about the outliers. Quantile regression cares about a specific slice of the distribution. The architecture stays fixed; the objective moves.

We work with the California Housing dataset throughout, 20,640 rows and 8 predictive features with MedHouseVal as our target in units of $100,000. It is a continuous regression problem with a heavy right tail and a hard cap near 5.0, which makes it ideal for exposing the differences between loss functions.

The data

Before any modeling, we inspect the raw data. The exploratory data analysis (EDA) reveals four findings that will shape everything that follows. First, MedHouseVal has a long right tail with 1,071 outliers by the 1.5 IQR rule, and the distribution is capped near 5.0. Second, MedInc has the strongest linear correlation with the target, which we will see dominate the coefficients in every model. Third, there are no missing values, so no imputation is needed. Fourth, AveOccup, Population, and AveRooms contain extreme outliers, which means scaling before any L1 or L2 penalized model is mandatory.

Figure 1 shows the target distribution.

The histogram makes the story concrete. The mass of districts sits between 1 and 3, but a long tail stretches toward the cap. A model that minimizes squared error will spend its capacity chasing those expensive districts, and that is exactly what we are about to see.

Linear regression family

We start with a naive baseline: predict the training mean for every test point. That gives a root mean squared error (RMSE) of 1.1449. Every model we build from here must beat that number, and the gap between our models and this baseline is the signal we are extracting from the features.

The first real model is Ordinary Least Squares, which minimizes squared error and has a closed-form solution through the normal equations. We solve it two ways, once with np.linalg.lstsq on the augmented design matrix and once with LinearRegression from scikit-learn.

# The normal equations solve least squares in closed form
X_b = np.hstack([np.ones((X_train_s.shape[0], 1)), X_train_s])
theta, _, _, _ = np.linalg.lstsq(X_b, y_train, rcond=None)
y_pred_ne = np.hstack([np.ones((X_test_s.shape[0], 1)), X_test_s]) @ theta

Both land at 0.7456, and the coefficients agree to within 1e-6. The closed form is elegant, but it does not scale to high dimensions or protect against correlated features. That is where Ridge regression enters. Ridge adds an L2 regularization penalty, shrinking coefficients toward zero without ever setting them exactly to zero. With alpha=1.0 we get 0.7456, essentially identical to OLS because the features are already standardized and the problem is well-conditioned.

Lasso regression takes a different route. Its L1 regularization and sparsity penalty drives some coefficients exactly to zero, which makes it a feature selector. With alpha=0.01 we get 0.7404, the best of the linear models so far, and exactly one coefficient is zeroed out. Elastic Net mixes both penalties with l1_ratio=0.5 and lands at 0.7416, between the two.