Linear and Logistic Regression

Linear and logistic regression are the workhorses of applied machine learning. They are the first models you should fit on almost any tabular problem, they are still used directly in credit scoring, medicine, advertising and economics because their coefficients can be read and audited, and their machinery (a weighted sum, a loss function, gradient descent, regularisation, softmax) is exactly what sits inside every neural network. If you understand this lesson, the output layer of a transformer holds no mysteries.

Intuition: the best straight edge

Imagine a scatter of points on graph paper and a ruler. Linear regression slides and tilts the ruler until the total squared vertical distance from the points to the ruler is as small as possible. Squaring makes large misses count much more than small ones, and it gives a smooth bowl-shaped error surface with one lowest point that we can find exactly.

Logistic regression uses the same ruler for a yes/no question. Instead of predicting the value directly, it computes a score along the ruler and squashes that score into a probability between 0 and 1. The ruler's line (where the probability is exactly 0.5) becomes a decision boundary separating the two classes.

Linear regression

The model

We have $n$ examples. Example $i$ has a feature vector $x_i \in \mathbb{R}^d$ and a real-valued target $y_i$. The model predicts

with weights (coefficients) $w \in \mathbb{R}^d$ and intercept (bias) $b$. It is convenient to append a constant 1 to every $x_i$ and fold $b$ into the weight vector; stacking all examples as rows gives a design matrix $X \in \mathbb{R}^{n \times (d+1)}$ and the predictions $\hat y = X\theta$, where $\theta = (b, w_1, \dots, w_d)$.

Least squares

We choose $\theta$ to minimise the mean squared error (MSE):

Squared error is not arbitrary: if the targets are generated as $y_i = \theta^\top x_i + \varepsilon_i$ with independent Gaussian noise $\varepsilon_i \sim \mathcal N(0, \sigma^2)$, minimising MSE is the same as maximising the likelihood of the data.

The normal equation

$J$ is a convex quadratic in $\theta$, so its minimum is where the gradient is zero:

These are the normal equations. If $X^\top X$ is invertible (no feature is an exact linear combination of others),

In code, never form the inverse explicitly; call a least-squares solver (np.linalg.lstsq), which uses numerically stable QR or SVD decompositions. Solving costs roughly $O(nd^2 + d^3)$, which is fine for thousands of features but not millions.

Gradient descent

For very large $n$ or $d$, or when no closed form exists (as in logistic regression), we use gradient descent: start from some $\theta$ and repeatedly step downhill,

where the learning rate $\eta > 0$ sets the step size. Too large and the iterates diverge; too small and progress crawls. Stochastic gradient descent (SGD) estimates the gradient from a small random batch of examples, trading noise for much cheaper steps. Gradient descent converges much faster when features are on similar scales, which is one reason to standardise them.

Worked example: three points

Data: $(x, y) = (1, 2), (2, 3), (3, 5)$. With one feature, the least-squares solution has the familiar form

• Means: $\bar x = 2$, $\bar y = 10/3 \approx 3.333$.
• Numerator: $(-1)(2 - 3.333) + (0)(3 - 3.333) + (1)(5 - 3.333) = 1.333 + 0 + 1.667 = 3$.
• Denominator: $1 + 0 + 1 = 2$.
• So $w = 1.5$ and $b = 3.333 - 1.5 \times 2 = 0.333$.

Predictions are $1.833, 3.333, 4.833$, residuals$0.167, -0.333, 0.167$.

One step of gradient descent from $w = b = 0$ with $\eta = 0.1$, for comparison: the residuals $\hat y - y$ are $(-2, -3, -5)$, so $\partial J/\partial w = \frac{2}{3}\sum (\hat y_i - y_i)x_i = \frac{2}{3}(-2 - 6 - 15) = -15.33$ and $\partial J/\partial b = \frac{2}{3}(-10) = -6.67$. The update gives $w = 1.533$, $b = 0.667$: already close to the answer, and further steps converge to $(1.5, 0.333)$.

Measuring fit: $R^2$

The coefficient of determination compares the model's squared error with that of always predicting the mean:

In the example, $\text{SSE} = 0.0278 + 0.1111 + 0.0278 = 0.1667$ and $\text{SST} = 1.778 + 0.111 + 2.778 = 4.667$, so $R^2 = 1 - 0.1667/4.667 = 0.964$: the line explains about 96% of the variance in $y$. On training data $R^2$ never decreases when you add features, so judge it on held-out data, where it can even be negative (worse than predicting the mean).

Assumptions, and which ones matter

Linear regression fits any data. The assumptions matter for whether the fit is good and whether standard errors and p-values are trustworthy:

1. Linearity: the expected target is linear in the features (after any transformations you apply). Check residual-vs-fitted plots for curves.
2. Independent errors: violated by time series and repeated measurements; standard errors become too small.
3. Constant error variance (homoscedasticity): violated when errors grow with the prediction, common for prices and counts; a log transform of $y$ often helps.
4. Normal errors: needed only for exact small-sample inference; irrelevant for prediction.
5. No perfect multicollinearity: if features are exact linear combinations, $X^\top X$ is singular. Strong (not perfect) correlation makes individual coefficients unstable even when predictions are fine.

For pure prediction, only (1) truly matters, and it is checked by validation error. For explanation, all of them do.

Interpreting coefficients

$w_j$ is the expected change in $\hat y$ for a one-unit increase in feature $j$, holding the other features fixed. Three cautions:

• Units matter. A coefficient of 0.001 on income in dollars is not "small". Standardise features (subtract the mean, divide by the standard deviation) to compare magnitudes; then $w_j$ is the change per standard deviation.
• "Holding others fixed" may be impossible in reality when features are correlated (you cannot raise square footage while holding room count fixed), and correlated features can split or flip each other's coefficients.
• Coefficients describe association, not causation. A positive coefficient on "number of firefighters" for fire damage does not mean firefighters cause damage.

Logistic regression

From scores to probabilities

For binary labels $y_i \in \{0, 1\}$ we want $P(y = 1 \mid x)$. A linear score $z = w^\top x + b$ can be any real number, so we pass it through the sigmoid (logistic) function:

$\sigma$ maps $(-\infty, \infty)$ to $(0, 1)$, with $\sigma(0) = 0.5$ and $\sigma(-z) = 1 - \sigma(z)$.

Odds and log-odds

The odds of an event with probability $p$ are $p/(1-p)$; a probability of 0.8 is odds of 4 (to 1). Inverting the sigmoid shows what logistic regression really models:

The log-odds (logit) are linear in the features. So a one-unit increase in $x_j$ adds $w_j$ to the log-odds, which multiplies the odds by $e^{w_j}$ (the odds ratio). With $w_j = 0.7$, each unit multiplies the odds by $e^{0.7} \approx 2.01$. Note that this is not a fixed change in probability: the same odds ratio moves $p$ from 0.5 to 0.67, but from 0.95 to only 0.97.

The decision boundary

Predicting class 1 when $\hat p \ge 0.5$ is the same as predicting it when $w^\top x + b \ge 0$. The boundary $w^\top x + b = 0$ is a straight line (a hyperplane in $d$ dimensions), so logistic regression is a linear classifier. Curved boundaries require transformed features (squares, interactions, splines). The 0.5 threshold is a choice: when false negatives cost more than false positives, lower it.

Log-loss

The squared error would make the optimisation non-convex here. Instead we maximise the likelihood $\prod_i \hat p_i^{\,y_i}(1-\hat p_i)^{1-y_i}$, which is equivalent to minimising the log-loss (binary cross-entropy):

It punishes confident mistakes severely: predicting $\hat p = 0.01$ for a true positive costs $-\log 0.01 = 4.6$, while $\hat p = 0.5$ costs only $0.69$. $L$ is convex, and its gradient has a strikingly simple form,

"prediction minus target, times input", the same shape as for linear regression. There is no closed-form solution, so we use iterative solvers (gradient descent, Newton's method, L-BFGS). If the classes are perfectly separable, the unregularised weights grow without bound; regularisation fixes this.

Worked example

A model has $w = 2$, $b = -3$ for a single feature. For $x = 2$: $z = 2(2) - 3 = 1$, so $\hat p = \sigma(1) = 1/(1 + e^{-1}) = 0.731$.

• If the true label is $y = 1$, the log-loss for this example is $-\log 0.731 = 0.313$.
• If $y = 0$, it is $-\log(1 - 0.731) = -\log 0.269 = 1.313$.
• The odds are $0.731/0.269 = e^{1} = 2.72$. Increasing $x$ by 1 multiplies the odds by $e^{2} = 7.39$, to $20.1$, i.e. $\hat p = 20.1/21.1 = 0.953$ (check: $\sigma(3) = 0.953$).
• The decision boundary is at $2x - 3 = 0$, i.e. $x = 1.5$.

Regularisation: ridge, lasso, elastic net

With many features, few examples or correlated features, least-squares and maximum-likelihood weights become large and unstable, and the model overfits. Regularisation adds a penalty on the size of the weights. For a data loss $\text{Loss}(w)$ (MSE or log-loss) and strength $\lambda \ge 0$:

The intercept is not penalised. (Exact scaling constants differ between textbooks and libraries; scikit-learn's Ridge minimises $\lVert y - Xw\rVert^2 + \alpha\lVert w\rVert^2$ without dividing by $n$, while its Lasso divides the squared error by $2n$.)

• Ridge shrinks all weights smoothly toward zero and has a closed form, $\hat w = (X^\top X + \lambda I)^{-1}X^\top y$ on centred data. Adding $\lambda I$ makes the matrix invertible even with collinear features. Among correlated features, ridge tends to share weight between them.
• Lasso drives some weights exactly to zero, performing feature selection. Geometrically, the L1 constraint region is a diamond whose corners lie on the axes, and the optimum often lands on a corner. Among a group of highly correlated features, lasso tends to pick one somewhat arbitrarily.
• Elastic net mixes both: sparse like lasso, but stable with correlated groups like ridge.

Regularisation penalises all weights equally, so features must be standardised first, or the penalty will hit small-unit features harder. Choose $\lambda$ by cross-validation. In scikit-learn, linear regressors use alpha for $\lambda$, but LogisticRegression uses C, the inverse strength: small C means strong regularisation. scikit-learn's LogisticRegression applies L2 regularisation with C=1.0 by default, which surprises people expecting a plain maximum-likelihood fit.

Multiclass classification

For $K > 2$ classes there are two standard approaches.

One-vs-rest (OvR) trains $K$ binary classifiers, class $k$ against all others, and predicts the class with the highest score. Simple and parallel, but the $K$ probabilities are not jointly normalised.

Multinomial (softmax) regression gives each class its own weight vector $w_k$ and bias $b_k$, computes scores $z_k = w_k^\top x + b_k$, and normalises them with the softmax:

Training minimises the cross-entropy $-\frac{1}{n}\sum_i \log P(y = y_i \mid x_i)$, which reduces to the log-loss when $K = 2$. For example, scores $z = (2.0, 1.0, 0.1)$ give $e^{z} = (7.39, 2.72, 1.11)$, summing to 11.21, so the probabilities are $(0.659, 0.242, 0.099)$; if the true class is the first, the loss is $-\log 0.659 = 0.417$. This exact computation is the final layer of nearly every neural network classifier, including the next-token prediction of language models. Current versions of scikit-learn's LogisticRegression use the multinomial formulation for multiclass problems with its default solver; wrap a model in OneVsRestClassifier if you specifically want OvR.

Code

First, three ways to fit the same linear regression, which should agree:

import numpy as np
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(0)
n = 200
x = rng.uniform(0, 10, size=(n, 1))
y = 3.0 + 2.0 * x[:, 0] + rng.normal(0, 1.0, size=n)   # true b=3, w=2

# Normal equation (use lstsq/solve, never an explicit inverse)
X = np.hstack([np.ones((n, 1)), x])                       # add bias column
theta_ne = np.linalg.lstsq(X, y, rcond=None)[0]

# Batch gradient descent on MSE
theta = np.zeros(2)
lr = 0.01
for _ in range(5000):
    grad = (2 / n) * X.T @ (X @ theta - y)
    theta -= lr * grad

sk = LinearRegression().fit(x, y)
print("normal eq :", theta_ne.round(3))
print("grad desc :", theta.round(3))
print("sklearn   :", round(sk.intercept_, 3), sk.coef_.round(3))
print("R^2       :", round(sk.score(x, y), 3))

All three recover $b \approx 3.05$, $w \approx 1.98$ (close to the true 3 and 2; the difference is noise), with $R^2 \approx 0.97$.

Next, regularised regression on scikit-learn's diabetes dataset (442 patients, 10 features, target is a disease progression score):

import numpy as np
from sklearn.datasets import load_diabetes
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet

X, y = load_diabetes(return_X_y=True)
models = {
    "OLS":        LinearRegression(),
    "Ridge a=10": Ridge(alpha=10),
    "Lasso a=1":  Lasso(alpha=1.0),
    "ENet a=1":   ElasticNet(alpha=1.0, l1_ratio=0.5),
}
for name, m in models.items():
    pipe = make_pipeline(StandardScaler(), m)
    r2 = cross_val_score(pipe, X, y, cv=5, scoring="r2").mean()
    coef = pipe.fit(X, y)[-1].coef_
    print(f"{name:11s} CV R2={r2:.3f}  zero coefs={np.sum(coef == 0)}/10")

Typical output: all of OLS, ridge and lasso reach a cross-validated $R^2$ of about 0.48, but lasso does it with 3 of the 10 coefficients set exactly to zero, a simpler model at no cost. The elastic net at this alpha is over-regularised ($R^2 \approx 0.46$): the strength must be tuned (for example with LassoCV, RidgeCV or ElasticNetCV), not guessed.

Finally, logistic regression with interpretable coefficients:

import numpy as np
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import log_loss, accuracy_score

data = load_breast_cancer()
X_tr, X_te, y_tr, y_te = train_test_split(
    data.data, data.target, test_size=0.25, stratify=data.target, random_state=0)

clf = make_pipeline(StandardScaler(), LogisticRegression(C=1.0, max_iter=1000))
clf.fit(X_tr, y_tr)
p = clf.predict_proba(X_te)[:, 1]
print("test log-loss: %.3f  accuracy: %.3f" % (log_loss(y_te, p), accuracy_score(y_te, p > 0.5)))

coef = clf[-1].coef_[0]
top = np.argsort(np.abs(coef))[::-1][:3]
for i in top:   # odds ratio per +1 standard deviation of the feature
    print(f"{data.feature_names[i]:22s} w={coef[i]:+.2f}  odds x{np.exp(coef[i]):.2f}")

This reaches about 0.96 test accuracy with a log-loss near 0.09. Because features were standardised, each coefficient is per standard deviation. In this dataset class 1 is benign, so the negative weights (for example on "radius error", about $-1.2$, odds ratio about 0.30) mean larger values make benign less likely: each extra standard deviation cuts the odds of benign to roughly a third. Remember the caveat from earlier: many of these 30 features are strongly correlated (radius, perimeter and area measure nearly the same thing), so individual coefficients should be interpreted cautiously.

Common pitfalls

• Not scaling before regularising, so the penalty depends on arbitrary units.
• Forgetting that LogisticRegression is regularised by default and that C is an inverse strength.
• Reading coefficients causally, or comparing raw coefficients across features with different units.
• Interpreting coefficients of strongly correlated features individually; they can split, shrink or flip sign.
• Reporting training $R^2$ as evidence of predictive power.
• Using linear regression for a binary target, which yields "probabilities" below 0 or above 1.
• Treating 0.5 as the only threshold; choose the threshold from the costs of errors.
• Computing the inverse $(X^\top X)^{-1}$ explicitly instead of using a stable solver.

In practice

• Fit a regularised linear or logistic model as your first serious baseline on every tabular problem. It trains in seconds, rarely fails badly and tells you whether the signal is roughly linear.
• Standardise numeric features and one-hot encode categoricals inside a pipeline; tune alpha or C on a logarithmic grid by cross-validation.
• Use L1 or elastic net when you want a sparse, explainable model; ridge when you want stability with many correlated features.
• For heavily skewed targets (prices, durations), model $\log y$; for counts, consider Poisson regression (PoissonRegressor).
• Logistic regression's probabilities are usually reasonably calibrated, which is why it remains standard in credit scoring and clinical risk scores; check calibration anyway (covered in model-evaluation-and-metrics).
• When the linear baseline is clearly beaten by tree ensembles (next lesson), that tells you the signal involves interactions and non-linearities, which you can either engineer as features or leave to a more flexible model.