What you’ll learn

• the difference between a cost function, a loss function and a metric
• the definition of MSE, and why we train on it but report RMSE
• how squaring changes what the model considers important
• MAE and Huber — what they optimise instead, and when each wins
• why the MSE surface is a convex bowl, and what that guarantees
• the statistical reason MSE is the default: it is maximum likelihood under Gaussian noise

Intuition

A model needs a single number that says how wrong it currently is — one number, because optimisation can only push one thing downhill. That number is the cost.

Choosing the cost is not a technicality. It is a statement of values. Squaring the errors declares that one miss of 10 is worse than ten misses of 1 (100 against 10), so the model will contort itself to avoid big misses even at the cost of many small ones. Taking absolute values declares they are equally bad. Both are defensible; they simply produce different models from identical data.

  flowchart LR
  X["Features X"] --> H["Model h(x)"]
  H --> P["Predictions ŷ"]
  Y["True targets y"] --> E["Errors ŷ - y"]
  P --> E
  E --> S["Square each one"]
  S --> M["Average → MSE"]
  M -->|"gradient"| H

Three words that get confused

The cost and the metric are allowed to differ, and often should. You train on MSE because it has clean gradients; you report RMSE because it is in the units of the target.

The math

Definition

Units are the square of the target’s units — squared dollars, squared degrees. Nobody has intuition for squared dollars, which is why we take the root:

Why train on MSE and report RMSE

The square root is monotonic, so whichever $\boldsymbol{\theta}$ minimises MSE also minimises RMSE. The minimiser is identical. But the gradient is not:

The RMSE gradient carries a $1/\sqrt{\text{MSE}}$ factor that blows up as the cost approaches zero. MSE’s gradient is clean and linear in the residual. Same answer, better-behaved arithmetic — so train on MSE, report RMSE.

Why the surface is a bowl

$\text{MSE}(\boldsymbol{\theta})$ is a quadratic form in $\boldsymbol{\theta}$ with Hessian

$\mathbf{X}^\top\mathbf{X}$ is positive semi-definite for any $\mathbf{X}$ — for any vector $\mathbf{v}$, $\mathbf{v}^\top\mathbf{X}^\top\mathbf{X}\mathbf{v} = \lVert\mathbf{X}\mathbf{v}\rVert^2 \geq 0$. So the cost is convex: no local minima, no saddle points, no bad starting position. Whatever optimiser you point at it will find the global minimum, and if the columns are linearly independent that minimum is unique.

figureMSE over the two parameters of a simple linear modelmatplotlib

The statistical justification

MSE is not an arbitrary choice. Assume the target is a linear function plus Gaussian noise:

The log-likelihood of the data is then

Only the last term depends on $\boldsymbol{\theta}$, and it is $-\frac{1}{2\sigma^2}$ times the sum of squared errors. Maximising the likelihood is exactly minimising MSE. Choosing squared error is the same act as assuming your noise is Gaussian — which also explains why MSE handles outliers badly: Gaussians have thin tails and treat a large deviation as nearly impossible.

Worked example by hand

The Phase 3 five-point set with the fitted line $\hat{y} = 0.6x + 2.2$:

Note that RMSE (0.693) exceeds MAE (0.64). That is not a coincidence: RMSE is always at least as large as MAE, and the gap widens as the errors become less uniform. A big RMSE-to-MAE ratio is a quick signal that a few large errors dominate.

Now break one point

Change the last observation from 5 to 25 and refit. The least-squares line becomes $\hat{y} = 4.6x - 5.8$:

MSE has jumped from 0.48 to 30.88, and — far worse — every single point is now badly predicted. The line abandoned four good observations to chase one bad one. The MAE-optimal line for the same data is $\hat{y} \approx 1.51x + 0.49$, which stays close to the original four points and simply accepts a large error on the outlier.

Comparing the three losses

Huber is quadratic near zero (so it has a smooth gradient where it matters) and linear in the tails (so one wild point cannot dominate). It is the pragmatic default when you suspect outliers but still want gradient-friendly behaviour.

figurePenalty as a function of a single residualmatplotlib

figureDragging one point upward, and watching each loss reactmatplotlib

Reading the plots

• The MAE line is flat. Not merely flatter — flat. Absolute-error fitting depends on the rank of the residuals, not their size, so once a point is above the line it does not matter whether it is above by 5 or by 500.
• The MSE line rises without limit. Each unit of extra distance adds a proportional pull on the slope. There is no point at which MSE stops caring.
• Huber sits between them, by construction. Below $\delta$ it behaves like MSE, above it like MAE, and $\delta$ is the dial.

See it move

Each pink square’s area is one point’s squared error. Watch the total area collapse as the line improves — that shrinking area is literally the cost.

From scratch

All three losses, plus the invariant that RMSE never falls below MAE:

import numpy as np
 
 
def mse(y, pred):
    return float(((pred - y) ** 2).mean())
 
 
def rmse(y, pred):
    return mse(y, pred) ** 0.5
 
 
def mae(y, pred):
    return float(np.abs(pred - y).mean())
 
 
def huber(y, pred, delta=1.0):
    e = np.abs(pred - y)
    quadratic = 0.5 * e**2
    linear = delta * (e - 0.5 * delta)
    return float(np.where(e <= delta, quadratic, linear).mean())
 
 
x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
y = np.array([2.0, 4.0, 5.0, 4.0, 5.0])
pred = 0.6 * x + 2.2
 
print(f"MSE   = {mse(y, pred):.4f}")     # MSE   = 0.4800
print(f"RMSE  = {rmse(y, pred):.4f}")    # RMSE  = 0.6928
print(f"MAE   = {mae(y, pred):.4f}")     # MAE   = 0.6400
print(f"Huber = {huber(y, pred):.4f}")   # Huber = 0.2400
 
print(rmse(y, pred) >= mae(y, pred))     # True — always

import numpy as np
 
 
def mse(y, pred):
    return float(((pred - y) ** 2).mean())
 
 
def rmse(y, pred):
    return mse(y, pred) ** 0.5
 
 
def mae(y, pred):
    return float(np.abs(pred - y).mean())
 
 
def huber(y, pred, delta=1.0):
    e = np.abs(pred - y)
    quadratic = 0.5 * e**2
    linear = delta * (e - 0.5 * delta)
    return float(np.where(e <= delta, quadratic, linear).mean())
 
 
x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
y = np.array([2.0, 4.0, 5.0, 4.0, 5.0])
pred = 0.6 * x + 2.2
 
print(f"MSE   = {mse(y, pred):.4f}")     # MSE   = 0.4800
print(f"RMSE  = {rmse(y, pred):.4f}")    # RMSE  = 0.6928
print(f"MAE   = {mae(y, pred):.4f}")     # MAE   = 0.6400
print(f"Huber = {huber(y, pred):.4f}")   # Huber = 0.2400
 
print(rmse(y, pred) >= mae(y, pred))     # True — always

With scikit-learn

import numpy as np
from sklearn.linear_model import LinearRegression, HuberRegressor
from sklearn.metrics import mean_squared_error, mean_absolute_error
 
X = np.array([1.0, 2.0, 3.0, 4.0, 5.0]).reshape(-1, 1)
clean = np.array([2.0, 4.0, 5.0, 4.0, 5.0])
dirty = np.array([2.0, 4.0, 5.0, 4.0, 25.0])   # last point corrupted
 
ols_clean = LinearRegression().fit(X, clean)
ols_dirty = LinearRegression().fit(X, dirty)
huber_dirty = HuberRegressor(epsilon=1.35).fit(X, dirty)
 
print(f"OLS on clean data : slope {ols_clean.coef_[0]:.2f}")    # slope 0.60
print(f"OLS on dirty data : slope {ols_dirty.coef_[0]:.2f}")    # slope 4.60
print(f"Huber on dirty    : slope {huber_dirty.coef_[0]:.2f}")  # slope 1.91
 
# Same predictions, two different metrics
pred = ols_dirty.predict(X)
print(f"MSE = {mean_squared_error(dirty, pred):.2f}")           # MSE = 30.88
print(f"MAE = {mean_absolute_error(dirty, pred):.2f}")          # MAE = 4.64

import numpy as np
from sklearn.linear_model import LinearRegression, HuberRegressor
from sklearn.metrics import mean_squared_error, mean_absolute_error
 
X = np.array([1.0, 2.0, 3.0, 4.0, 5.0]).reshape(-1, 1)
clean = np.array([2.0, 4.0, 5.0, 4.0, 5.0])
dirty = np.array([2.0, 4.0, 5.0, 4.0, 25.0])   # last point corrupted
 
ols_clean = LinearRegression().fit(X, clean)
ols_dirty = LinearRegression().fit(X, dirty)
huber_dirty = HuberRegressor(epsilon=1.35).fit(X, dirty)
 
print(f"OLS on clean data : slope {ols_clean.coef_[0]:.2f}")    # slope 0.60
print(f"OLS on dirty data : slope {ols_dirty.coef_[0]:.2f}")    # slope 4.60
print(f"Huber on dirty    : slope {huber_dirty.coef_[0]:.2f}")  # slope 1.91
 
# Same predictions, two different metrics
pred = ols_dirty.predict(X)
print(f"MSE = {mean_squared_error(dirty, pred):.2f}")           # MSE = 30.88
print(f"MAE = {mean_absolute_error(dirty, pred):.2f}")          # MAE = 4.64

One corrupted observation multiplies the OLS slope by more than seven, from 0.60 to 4.60. Huber still drifts — 1.91 — because five points is a brutal test and one of them is now an extreme value, but it absorbs less than half the damage OLS does.

Pitfalls

🧪 Try It Yourself

Exercise 1 – Compute MSE by hand

Exercise 2 – RMSE is never below MAE

Exercise 3 – Implement the Huber loss

Exercise 4 – Mean minimises MSE, median minimises MAE

Exercise 5 – Huber resists what OLS cannot

Recap

• Loss is per instance, cost is over the training set, and the metric is what you report.
• $\text{MSE} = \frac{1}{m}\sum e_i^2$; train on it, report RMSE because it is in target units.
• The cost surface is convex because $\mathbf{X}^\top\mathbf{X}$ is positive semi-definite — one global minimum, no traps.
• Minimising MSE is maximum likelihood under Gaussian noise, which is also why it fears outliers.
• MSE targets the mean, MAE the median, Huber the mean with a cap on how much any one point can pull.
• On the five-point set: MSE 0.48, RMSE 0.693, MAE 0.64, Huber 0.24. Corrupt one point and MSE jumps to 30.88 while the Huber slope barely moves.

Exercise 6 – Which statistic is each loss estimating?

Next

Continue to Gradient Descent Explained — how to actually walk down the bowl when the closed-form solution is too expensive to compute.