Linear Regression: Predicting Outcomes

Linear regression is the first algorithm most people learn and the one many teams still ship. It predicts a number from other numbers by adding up weighted inputs. That sounds too simple to be useful, yet it forecasts demand, prices houses, estimates delivery times and sits underneath plenty of dashboards labelled “AI”. It trains in milliseconds, explains itself through its weights and gives you a baseline that any fancier model has to beat.

This article covers what the model actually is, how it finds its weights, how to train and evaluate one in Python, which assumptions matter in practice and what to do when they break. It assumes you know the vocabulary from key concepts and terminologies in ML: features, labels, training and test sets, loss functions.

What linear regression predicts

Linear regression is a supervised learning method for regression problems, where the label is a continuous number rather than a category. Given examples with known answers, it learns a formula that maps inputs to an output, then applies that formula to inputs it has never seen.

The running example in this article is a monthly cloud bill. Each row is one month for one team, with three features:

• vCPU-hours of compute used
• average gigabytes of storage held
• gigabytes of data sent out to the internet (egress)

The label is the bill in dollars. The data is synthetic and the rates are made up, which is the point: because you know the true formula that generated it, you can check whether the model recovers it.

The model in one line

With a single feature, linear regression fits a straight line:

ŷ = w · x + b

Here x is the input, w is the weight (the slope), b is the intercept (the value when x is zero) and ŷ is the prediction. With several features it becomes a weighted sum:

ŷ = w₁·x₁ + w₂·x₂ + w₃·x₃ + b

For the cloud bill, w₁ is the cost per vCPU-hour, w₂ the cost per stored gigabyte, w₃ the cost per egress gigabyte and b any fixed monthly charge. The weights and intercept are the model’s parameters. Training means finding values for them that make the predictions close to the real bills.

flowchart LR
    X1["vCPU-hours"] -->|"× w1"| S["Sum + b"]
    X2["Storage GB"] -->|"× w2"| S
    X3["Egress GB"] -->|"× w3"| S
    S --> Y["Predicted bill"]

“Linear” refers to the parameters, not the shape of the curve. A model like ŷ = w₁·x + w₂·x² + b is still linear regression, because the prediction is a weighted sum of inputs; one of the inputs just happens to be x². That detail is what lets linear regression fit curves, as you will see later.

How the model learns its weights

For any choice of weights, each training example has a residual: the gap between the real value and the prediction, y − ŷ. Some residuals are positive, some negative. Linear regression picks the weights that minimize the sum of squared residuals, which is why the standard method is called ordinary least squares (OLS).

Dividing that sum by the number of examples gives the mean squared error (MSE), the loss function you will see in most libraries:

MSE = (1/n) · Σ (yᵢ − ŷᵢ)²

Squaring does two jobs. It stops positive and negative errors from cancelling out, and it punishes large errors much more than small ones: a miss of 10 costs 100, a miss of 2 costs only 4. That second property is a strength when big misses are what you care about and a weakness when your data has outliers, because a single extreme point can drag the whole line towards it.

There are two common ways to find the weights that minimize MSE.

The closed-form solution

Least squares has an exact answer. If you put the features in a matrix X (with a column of ones for the intercept) and the labels in a vector y, the best weights solve the normal equation:

w = (Xᵀ X)⁻¹ Xᵀ y

You should not compute that inverse directly in code, because it is numerically fragile when features are correlated. Libraries use more stable decompositions instead. NumPy’s lstsq solves the same problem safely. Save the data generator below as data.py; the later examples import from it:

import numpy as np

rng = np.random.default_rng(42)
n = 500

vcpu_hours = rng.uniform(200, 5000, n)     # compute used in the month
storage_gb = rng.uniform(50, 2000, n)      # average storage held
egress_gb = rng.uniform(0, 1500, n)        # data sent out to the internet

# Made-up rates plus a fixed monthly charge and some noise
bill = (
    40
    + 0.045 * vcpu_hours
    + 0.02 * storage_gb
    + 0.09 * egress_gb
    + rng.normal(0, 15, n)
)

X = np.column_stack([vcpu_hours, storage_gb, egress_gb])
y = bill

Then fit it:

import numpy as np
from data import X, y

# Add a column of ones so the intercept is learned like any other weight
X_b = np.column_stack([np.ones(len(X)), X])

# Solve the least squares problem directly
w, *_ = np.linalg.lstsq(X_b, y, rcond=None)
print("intercept:", round(w[0], 2))
print("weights:  ", np.round(w[1:], 4))

intercept: 35.81
weights:   [0.0447 0.022  0.0903]

The true rates were 0.045, 0.02 and 0.09, with a fixed charge of 40. The model recovered the rates closely. The intercept is less precise because it is an extrapolation to zero usage, a point no month in the data comes near.

Gradient descent

The closed form needs the whole dataset in memory and gets expensive as the number of features grows. Gradient descent takes the other route: start with any weights, measure which direction reduces the loss, take a small step that way and repeat.

flowchart TD
    A["Start with zero weights"] --> B["Predict for all rows"]
    B --> C["Compute MSE and its gradient"]
    C --> D["Step weights against the gradient"]
    D --> E{"Loss still falling?"}
    E -->|"Yes"| B
    E -->|"No"| F["Stop: weights have converged"]

The gradient of MSE with respect to the weights has a simple form, so the whole loop fits in a few lines:

import numpy as np
from data import X, y

# Standardize features so one learning rate suits every weight
mean, std = X.mean(axis=0), X.std(axis=0)
X_s = (X - mean) / std

w = np.zeros(X_s.shape[1])
b = 0.0
lr = 0.1

for step in range(500):
    y_hat = X_s @ w + b
    error = y_hat - y
    grad_w = 2 * X_s.T @ error / len(y)
    grad_b = 2 * error.mean()
    w -= lr * grad_w
    b -= lr * grad_b

# Convert back to the original units to compare with lstsq
w_orig = w / std
b_orig = b - (w_orig * mean).sum()
print("intercept:", round(b_orig, 2))
print("weights:  ", np.round(w_orig, 4))

intercept: 35.81
weights:   [0.0447 0.022  0.0903]

Same answer, reached iteratively. The standardization step matters. vCPU-hours run into the thousands while the other features are smaller, so without scaling, a learning rate small enough for the large feature crawls on the small ones, and one large enough for the small ones diverges on the large one.

For most tabular problems that fit in memory, the closed form is the right choice and the one scikit-learn uses. Gradient descent matters because it is the same procedure that trains neural networks; linear regression is simply the smallest model you can watch it work on.

Training and evaluating with scikit-learn

In practice you will use a library. scikit-learn’s LinearRegression wraps a least squares solver and follows the same fit and predict interface as every other model in the library, which makes swapping models later trivial.

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_error, root_mean_squared_error, r2_score
from sklearn.model_selection import train_test_split
from data import X, y

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=0
)

model = LinearRegression()
model.fit(X_train, y_train)

pred = model.predict(X_test)
print("intercept:", round(model.intercept_, 2))
print("weights:  ", np.round(model.coef_, 4))
print("MAE: ", round(mean_absolute_error(y_test, pred), 2))
print("RMSE:", round(root_mean_squared_error(y_test, pred), 2))
print("R^2: ", round(r2_score(y_test, pred), 4))

# Predict one new month: 3000 vCPU-hours, 800 GB stored, 400 GB egress
print("forecast:", round(model.predict([[3000, 800, 400]])[0], 2))

intercept: 37.86
weights:   [0.0445 0.0219 0.0885]
MAE:  13.09
RMSE: 15.64
R^2:  0.9546
forecast: 224.38

The weights differ slightly from the earlier fit because the model now trains on 80 percent of the rows. The metrics are computed on the 20 percent it never saw, which is the only honest way to judge a model; the model training and evaluation article explains why.

If you are following an older tutorial, you may see load_boston used as the example dataset. It was removed in scikit-learn 1.2, and the library now points to the California housing dataset (fetch_california_housing) instead. root_mean_squared_error is also fairly recent; on older versions you take the square root of mean_squared_error yourself.

Reading the metrics

Each metric answers a different question, so report more than one.

RMSE close to the noise level you put in (a standard deviation of 15) is a good sign: the model has captured the signal, and what remains is randomness no model could predict. On real data you do not know the noise level, so compare against a baseline instead, such as always predicting last month’s bill.

A high R² does not mean the model is right. It can be high on a model that is badly wrong in the region you care about, and it always rises when you add features, even useless ones. Treat it as a summary, not a verdict.

Interpreting the weights

The weights are why people reach for linear regression when they need to explain a prediction. Each one says how much the prediction changes when that feature rises by one unit and the others stay fixed. Here every extra egress gigabyte adds about nine cents to the predicted bill.

Two cautions apply. First, the size of a weight depends on the feature’s units, so you cannot compare weights directly to rank importance unless the features are standardized. Second, “the others stay fixed” is doing a lot of work. When features move together, the model cannot tell them apart, which is the multicollinearity problem in the next section.

The assumptions that matter

Linear regression always produces an answer. Whether that answer can be trusted depends on a handful of assumptions. Some affect the predictions; others mainly affect the confidence intervals and p-values that statistics packages report.

The first check to make on any fitted model is a residual plot: predicted value on the x-axis, residual on the y-axis. A healthy model shows a shapeless cloud around zero. A curve means the relationship is not linear. A funnel means the error grows with the size of the prediction, which is common for money, where large bills vary by larger amounts.

flowchart TD
    A["Fit the model"] --> B["Plot residuals vs predictions"]
    B --> C{"What shape?"}
    C -->|"Random cloud"| D["Assumptions look fine"]
    C -->|"Curve"| E["Add terms or transform features"]
    C -->|"Funnel"| F["Log-transform the label"]
    E --> A
    F --> A

Multicollinearity in practice

Suppose you add a fourth feature, instance-hours, that is almost exactly vCPU-hours divided by four, because every instance in the fleet has four vCPUs. The two features carry the same information. Least squares can trade weight between them freely, so small changes in the training data produce very different weights:

import numpy as np
from sklearn.linear_model import LinearRegression
from data import X, y, rng

# A second compute feature that is almost exactly vCPU-hours / 4
instance_hours = X[:, 0] / 4 + rng.normal(0, 1, len(X))
X2 = np.column_stack([X, instance_hours])

for seed in [1, 2, 3]:
    idx = np.random.default_rng(seed).choice(len(X2), 300, replace=False)
    ols = LinearRegression().fit(X2[idx], y[idx])
    print(f"sample {seed}: vcpu={ols.coef_[0]:+.4f}  instance={ols.coef_[3]:+.4f}")

sample 1: vcpu=+0.1525  instance=-0.4324
sample 2: vcpu=+0.3380  instance=-1.1745
sample 3: vcpu=+0.0792  instance=-0.1366

The true cost per vCPU-hour is 0.045, yet the fitted weights swing wildly and instance-hours gets a negative weight, as if running instances made the bill cheaper. The predictions are still fine, because the two errors cancel. The interpretation is useless. If you only need predictions, this may not matter; if anyone will read the weights, drop one of the duplicates or regularize.

Beyond a straight line

Two extensions cover most of what plain linear regression cannot do on its own.

Polynomial and engineered features

Because “linear” refers to the weights, you can fit curves by feeding the model transformed inputs: squares, logs, interactions such as storage_gb × egress_gb, or flags for categories. scikit-learn’s PolynomialFeatures generates squares and interactions automatically. The trade-off is that the number of features grows quickly, and with it the risk of overfitting, where the model memorizes noise in the training set.

Regularization

Regularization adds a penalty on the size of the weights to the loss, so the model prefers smaller, more stable weights unless the data strongly supports larger ones. The strength of the penalty is a hyperparameter, called alpha in scikit-learn.

Penalties depend on the scale of the weights, so standardize features before regularizing. A pipeline keeps the scaler inside cross-validation, which avoids leaking test-set statistics into training:

from sklearn.linear_model import Ridge
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from data import X, y

for alpha in [0.1, 1.0, 10.0, 100.0]:
    model = make_pipeline(StandardScaler(), Ridge(alpha=alpha))
    scores = cross_val_score(model, X, y, cv=5, scoring="r2")
    print(f"alpha={alpha:>6}: mean R^2 = {scores.mean():.4f}")

alpha=   0.1: mean R^2 = 0.9560
alpha=   1.0: mean R^2 = 0.9560
alpha=  10.0: mean R^2 = 0.9553
alpha= 100.0: mean R^2 = 0.9130

On this clean dataset a small penalty changes nothing and a large one hurts, because there is nothing to fix. On real data with many correlated features, a moderate alpha often generalizes better than plain least squares. RidgeCV and LassoCV search over alpha for you.

Linear regression in the cloud

You do not need a notebook to train a linear model. The major platforms offer it close to where the data already lives.

BigQuery ML trains models with SQL, inside the warehouse. A linear regression is a CREATE MODEL statement with model_type = 'linear_reg':

CREATE OR REPLACE MODEL `my_dataset.bill_model`
OPTIONS (
  model_type = 'linear_reg',
  input_label_cols = ['bill']
) AS
SELECT vcpu_hours, storage_gb, egress_gb, bill
FROM `my_dataset.monthly_usage`;

SELECT * FROM ML.EVALUATE(MODEL `my_dataset.bill_model`);

SELECT * FROM ML.PREDICT(
  MODEL `my_dataset.bill_model`,
  (SELECT 3000 AS vcpu_hours, 800 AS storage_gb, 400 AS egress_gb)
);

Amazon SageMaker AI ships a built-in Linear Learner algorithm for regression and classification on large datasets, and Azure Machine Learning includes linear regression in its designer and automated ML. For a model this small, the deciding factor is usually where the data is and where predictions are needed, not training speed. If you are setting up an environment from scratch, the simple ML project article covers the local tooling.

When to use it, and when not to

Linear regression is the right first model when:

• the label is a continuous number and you need a fast, honest baseline
• stakeholders need to understand why a prediction came out the way it did
• data is limited, since a model with few parameters overfits less
• predictions must be cheap to compute or easy to embed, such as in SQL or a spreadsheet

Look elsewhere when the relationship is strongly non-linear with interactions you cannot engineer by hand, when the label is a category (that is logistic regression territory, the next article in this series), or when the data is images, audio or text. Even then, fit a linear model first. If a gradient-boosted ensemble only beats it by a small margin, the simpler model may be the better product.

For more depth, the scikit-learn linear models guide covers every variant mentioned here, and the BigQuery ML documentation covers training models in SQL.

Key takeaways

• Linear regression predicts a number as a weighted sum of features plus an intercept; “linear” refers to the weights, so transformed features let it fit curves.
• Ordinary least squares picks the weights that minimize squared error, either exactly with a solver or iteratively with gradient descent, which needs scaled features.
• Evaluate on held-out data with more than one metric: MAE for the typical miss, RMSE for large misses, R² as a summary.
• Weights are interpretable per unit of each feature, but only compare in size after standardization, and they become unstable when features are strongly correlated.
• Check a residual plot before trusting the model; curves and funnels point to specific fixes.
• Ridge, Lasso and ElasticNet add a penalty that stabilizes weights or selects features; standardize first and tune alpha with cross-validation.