MSE: the mean squared error

Why we need it

If a model predicts a number, it was almost certainly trained by reducing MSE, the mean squared error. It is the “default” loss for regression: it has an exact solution for linear models, a smooth gradient for neural networks and a clear statistical meaning. The method of least squares, published by Legendre in 1805 and tied to the normal law of errors by Gauss in 1809, is exactly the minimization of MSE.

Where it appears:

• linear regression — ordinary least squares (OLS) with a one-formula solution;
• neural networks for regression — MSELoss in PyTorch, "mse" in Keras; autoencoders that reconstruct their input;
• diffusion models — the network learns to predict the added noise, and the loss is the MSE between the true and the predicted noise;
• reinforcement learning — the value function is often trained on the squared TD error;
• gradient boosting — loss="squared_error": the negative gradient of MSE is the residual, so every tree is fitted to the residuals;
• signal and image processing — PSNR is computed from MSE, and the Kalman filter minimizes the expected squared error.

Intuition: the area of the squares

The word “squared” can be taken literally here. Draw a square with side $|e_i|$ on each error. Its area is $e_i^2$, and MSE is the average area of the squares. Compare with MAE, where lengths of segments are averaged: area grows as the square of the side, so big errors take up a disproportionate amount of space.

Below is the cheat sheet’s running example: $y=(3,\,5,\,8)$, $\hat y=(2.5,\,5,\,10)$. Drag the blue predictions. On the lower plot the errors sit on the parabola $e^2$; MSE is the average height of these points.

The formula, bone by bone

The same four steps as in MAE, only with a square instead of the absolute value.

$y_i$

the true value for object $i$;

$\hat y_i$

the model’s prediction;

$e_i=y_i-\hat y_i$

the error (residual); its sign disappears after squaring;

$e_i^2$

the square: non-negative, smooth, and growing faster the larger the error;

$\sum e_i^2$

the sum of squared residuals — SSE (sum of squared errors), or RSS (residual sum of squares);

$\frac1n$

averaging: MSE does not depend on the sample size.

Why a square and not the absolute value? There are three reasons. Smoothness: $e^2$ has the derivative $2e$ everywhere, including at zero. Algebra: a sum of squares is a quadratic form, and its minimum is found by solving a linear system (the normal equations). Geometry: $\sum e_i^2$ is a squared Euclidean distance, and the best prediction is an orthogonal projection, as in Pythagoras’ theorem.

In vector form MSE is the squared Euclidean norm of the error vector divided by $n$:

Hence its second name, the L2 loss. The units of MSE are the units of $y$ squared: roubles², minutes². To get back to the original units, take the square root and you get RMSE; in the example $\sqrt{1.417}\approx1.190$.

Properties and the best constant

• Non-negativity: $\mathrm{MSE}\ge0$, with equality only for exact predictions.
• Quadratic scaling: multiply $y$ and $\hat y$ by $k$ and MSE is multiplied by $k^2$. Converting roubles to thousands of roubles shrinks MSE a million times.
• Smoothness and convexity: MSE is a quadratic function of the predictions; for a linear model, a quadratic function of the parameters. If the columns of the feature matrix (together with the column of ones for the intercept) are linearly independent, it is strictly convex and the minimum is unique.
• Error decomposition: $\mathrm{MSE}=\bar e^{\,2}+\operatorname{Var}(e)$ — the squared mean bias plus the spread of the errors. In the example $\bar e=-0.5$: $1.417=0.25+1.167$.
• Comparison with MAE: $\mathrm{MSE}\ge\mathrm{MAE}^2$ (the mean of squares is at least the square of the mean). In the example: $1.417\ge0.833^2\approx0.694$.

The best constant is the mean

If the model answers with a single number $c$, MSE as a function of $c$ is a parabola. Expand it by adding and subtracting the mean $\bar y$:

The first term does not depend on $c$; the second is non-negative and vanishes only at $c=\bar y$. So the minimum of MSE is unique and is attained at the mean, and its value is the variance of the data. The derivative $2(c-\bar y)$ is a straight line proportional to the distance from the mean.

The derivative: a force proportional to the error

The derivative of the square is $2e$. The further the prediction is from the truth, the harder the gradient pulls it back, and at the minimum itself the force fades out smoothly. This is exactly what makes MSE so convenient for gradient descent.

Look at gradient descent on a single error. The step $\hat y\leftarrow\hat y-\eta\,\partial e^2/\partial\hat y=\hat y+2\eta e$ turns the error into $e\leftarrow(1-2\eta)\,e$. Each step multiplies the error by the same number — a geometric progression.

The bowl over the parameters: least squares and gradient descent

Take the line $\hat y=wx+b$ and five points: $x=(1,\,2,\,3,\,4,\,5)$, $y=(1.3,\,2,\,3.1,\,4.5,\,7)$. MSE as a function of $(w,\,b)$ is a quadratic form, and its graph is an elliptic bowl (a paraboloid) with a single bottom.

The exact solution: the normal equations

Setting the gradient to zero gives a linear system — the normal equations — and the least-squares solution:

For our points $w^*=1.39$, $b^*=-0.59$, $\mathrm{MSE}^*=0.2294$. In practice the inverse matrix is not computed: QR or SVD decompositions are used (np.linalg.lstsq), and on large data, gradient descent.

Gradient descent and the shape of the bowl

The speed of descent is set by the shape of the bowl — the eigenvalues of the Hessian $H=\frac2nX^\top X$. Along the steep direction ($\lambda_{\max}$) too large a step throws you over the bottom; along the flat one ($\lambda_{\min}$) the descent crawls. For our data $\lambda_{\max}=23.66$ and $\lambda_{\min}=0.338$; their ratio, the condition number, is about $70$: the bowl is stretched into a narrow ravine, and any step above $\eta_{\text{crit}}=2/\lambda_{\max}\approx0.0845$ diverges.

Bias and variance

MSE on new data splits into three parts, and this is perhaps the most important formula about overfitting. Let $y=f(x)+\varepsilon$, where the noise has mean $0$ and variance $\sigma^2$, and let the model $\hat f$ be trained on a random sample. Then the expected error at a point $x$ is:

$\mathbb E\hat f(x)$

the average prediction over all possible training samples;

$\text{bias}^2$

how far the average model misses the truth — the price of a too simple model;

$\text{variance}$

how much the model changes from sample to sample — the price of a too flexible model;

$\sigma^2$

irreducible noise: no model can remove it.

The probabilistic meaning: Gaussian noise

Why the square of all things? The classical answer was given by Gauss (1809). Let $y_i=\hat y_i+\varepsilon_i$, where the noise $\varepsilon_i$ is independent and normally distributed, $N(0,\sigma^2)$. The density of one observation is a bell around the prediction, and the likelihood of the sample is the product of the heights of the bells at the data points. The negative log-likelihood is:

For any fixed $\sigma$ this is an increasing linear function of MSE. Hence minimizing MSE is maximum likelihood under Gaussian noise. Minimizing over $\sigma$ gives one more link: $\hat\sigma^2=\mathrm{MSE}$. Later (1821–1823) Gauss gave a second justification that needs no normality: if the errors are unbiased, uncorrelated and have equal variance, the least-squares estimate has the smallest variance among all linear unbiased estimates (the Gauss–Markov theorem).

One more statement of the same fact: among all functions $g(x)$ the expected squared error $\mathbb E(y-g(x))^2$ is smallest for $g(x)=\mathbb E[y\mid x]$. A model trained on MSE estimates the conditional mean. If the noise is asymmetric or heavy-tailed, the conditional mean can describe the “typical” outcome poorly — then one looks at the median (MAE) or at quantiles.

Outliers: unbounded influence

The flip side of the square is sensitivity to anomalous points. The gradient from an object is proportional to its error, so a distant outlier pulls the model harder the further away it is. Robust statistics says that least squares has an unbounded influence function and a breakdown point of zero: a single point is enough to drive the line arbitrarily far.

What to do about outliers if you need MSE:

• replace the loss with Huber: a square for small errors and the absolute value for large ones — the smoothness of MSE and the bounded influence of MAE;
• transform the target ($\log y$ for prices and incomes) — but then the model estimates the mean of the logarithm, not the mean of $y$;
• find and inspect outliers by their residuals: points with a huge $e_i^2$ stand out immediately;
• use robust fitting algorithms (RANSAC, Theil–Sen) that learn from subsets of the data.

In code

MSE and its gradient in NumPy, the exact least-squares solution, and a check that the best constant is the mean:

pythonimport numpy as np

y     = np.array([3.0, 5.0, 8.0])
y_hat = np.array([2.5, 5.0, 10.0])

def mse(y, y_hat):
    return np.mean((y - y_hat) ** 2)

def mse_grad(y, y_hat):
    # gradient w.r.t. the predictions: 2(ŷ − y) / n
    return 2 * (y_hat - y) / y.size

print(mse(y, y_hat))        # 1.4166666666666667
print(mse_grad(y, y_hat))   # [-0.33333333  0.          1.33333333]

# the least-squares line: the exact solution of the normal equations
x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
t = np.array([1.3, 2.0, 3.1, 4.5, 7.0])
A = np.column_stack([x, np.ones_like(x)])
w, b = np.linalg.lstsq(A, t, rcond=None)[0]
print(round(w, 4), round(b, 4), round(mse(t, A @ [w, b]), 4))   # 1.39 -0.59 0.2294

# the best constant under MSE is the mean; MSE(c) = Var + (c − ȳ)²
data = np.array([2.0, 3.0, 5.0, 9.0, 16.0])
print(data.mean(), mse(data, data.mean()), mse(data, 12.0))   # 7.0 26.0 51.0

In PyTorch it is MSELoss. The gradient with respect to the predictions is $2(\hat y-y)/n$; descent with step $0.05$ on the raw feature from Figure 4 reaches the same least-squares solution, but slowly: to two decimals it takes about 300 steps (the listing runs 2000 to be safe):

pythonimport torch
import torch.nn as nn

y     = torch.tensor([3.0, 5.0, 8.0])
y_hat = torch.tensor([2.5, 5.0, 10.0], requires_grad=True)

loss = nn.MSELoss()(y_hat, y)    # same as ((y_hat - y) ** 2).mean()
loss.backward()
print(loss.item())               # 1.4166666269302368
print(y_hat.grad)                # tensor([-0.3333,  0.0000,  1.3333])

# gradient descent for a line: with η > 2/λ_max ≈ 0.0845 it diverges
x = torch.tensor([1.0, 2.0, 3.0, 4.0, 5.0])
t = torch.tensor([1.3, 2.0, 3.1, 4.5, 7.0])
w = torch.zeros(1, requires_grad=True)
b = torch.zeros(1, requires_grad=True)
opt = torch.optim.SGD([w, b], lr=0.05)
for step in range(2000):
    opt.zero_grad()
    nn.functional.mse_loss(w * x + b, t).backward()
    opt.step()
print(round(w.item(), 2), round(b.item(), 2))   # 1.39 -0.59

scikit-learn has MSE and RMSE metrics (the regressors’ score method returns $R^2$, though), LinearRegression minimizes MSE exactly, and boosting uses it as the default loss:

pythonimport numpy as np
from sklearn.metrics import mean_squared_error, root_mean_squared_error
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import HistGradientBoostingRegressor

print(mean_squared_error([3, 5, 8], [2.5, 5, 10]))        # 1.4166666666666667
print(root_mean_squared_error([3, 5, 8], [2.5, 5, 10]))   # 1.1902380714238083

X = np.array([[1.0], [2.0], [3.0], [4.0], [5.0]])
y = np.array([1.3, 2.0, 3.1, 4.5, 7.0])
lr = LinearRegression().fit(X, y)                    # minimizes MSE exactly
print(lr.coef_.round(4), round(lr.intercept_, 4))    # [1.39] -0.59

# boosting on MSE: each tree is fitted to the residuals y − F(x)
gbr = HistGradientBoostingRegressor(loss="squared_error")

Common mistakes

1. Reading MSE as the “average error”. MSE $=1.417$ in the example is in squared units; the average error is of order $\sqrt{1.417}\approx1.19$ (RMSE) or $0.833$ (MAE).
2. MSE on data with outliers. A few anomalies will drag the model along (Figure 7). Clean the data, or use Huber or MAE.
3. Unscaled features. A stretched bowl means slow descent and a small admissible step (Figure 4). Standardize the features.
4. Too large a step. With $\eta>2/\lambda_{\max}$ the loss grows exponentially and quickly turns into inf and nan.
5. MSE for classification on top of a sigmoid. The gradient $2(p-y)\,p(1-p)$ vanishes when the model is confidently wrong ($p\approx0$ for $y=1$). For probabilities use BCE or cross-entropy.
6. Judging a model by its training MSE. A low MSE on the training data can mean memorized noise (Figure 5, $d=10$). Look at validation.

Links to other formulas

• MAE: the absolute value instead of the square — the median instead of the mean, robustness to outliers, but a kink at zero.
• RMSE: $\sqrt{\mathrm{MSE}}$ — the same optimal model, but the error in the units of $y$.
• Huber: MSE near zero and MAE in the tails.
• BCE: the analogue of MSE for probabilities — also maximum likelihood, but for the Bernoulli distribution.
• SGD and Adam: MSE shows clearly how the step and the conditioning of the problem determine convergence.

Summary

Check yourself