Inventing a Loss Function

The sum is -0.5, which looks tiny. But look at house 2: the prediction is off by 0.9, and house 5 is off by 1.0. The positive and negative errors cancel each other out. We need a formula where errors can't cancel.

actual    = [1.0, 1.1, 0.9, 0.5, 3.0]
predicted = [0.9, 2.0, 0.7, 0.3, 2.0]

total_error = sum(p - a for a, p in zip(actual, predicted))
print(f"Sum of differences: {total_error}")
# Output: -0.5
# This is misleading — every prediction is wrong, but the sum is near zero.

Idea 1 — Absolute value: Take |predicted - actual| for each house. Negatives become positive, so nothing cancels.

Idea 2 — Squaring: Take (predicted - actual)² for each house. Any number squared is positive, so nothing cancels.

These are exactly the two approaches the ML community uses!

MSE jumps from 0.49 to 1882 — a 3800× increase. MAE goes from 0.56 to 20 — about 36×. MSE is far more sensitive to outliers because squaring amplifies large errors.

actual    = [1.0, 1.1, 0.9, 0.5, 3.0]
predicted = [0.9, 2.0, 0.7, 0.3, 2.0]

print("Normal data:")
print(f"  MSE: {compute_mse(actual, predicted):.4f}")   # 0.488
print(f"  MAE: {compute_mae(actual, predicted):.4f}")   # 0.560

predicted_outlier = [0.9, 2.0, 0.7, 0.3, 100.0]
print("\nWith outlier (house 5 predicted=100):")
print(f"  MSE: {compute_mse(actual, predicted_outlier):.4f}")   # 1882.088
print(f"  MAE: {compute_mae(actual, predicted_outlier):.4f}")   # 19.960

MSE: Squaring amplifies large errors disproportionately. One outlier can dominate the entire loss, pulling the model toward that single point.

MAE: The function |x| has a sharp corner (a "kink") at x = 0. This will turn out to be a problem when we try to compute gradients — which is essential for training ML models.

The gradient of $\text{diff}^2$ is $2 \cdot \text{diff}$. Notice that when the error is large (e.g. diff=3), the gradient is also large (6). The gradient grows without bound as the error grows — this is why MSE reacts so strongly to outliers.

def squared_gradient(diff):
    return 2 * diff

for d in [-2, -1, 0, 1, 3]:
    analytical = squared_gradient(d)
    numerical  = numerical_gradient(squared_loss, d)
    print(f"diff={d:+}: analytical={analytical:+.4f}, numerical={numerical:+.6f}")

The gradient is +1 for positive diff, -1 for negative diff. Unlike the squared loss, the gradient doesn't grow with the error — it's always ±1 regardless of how big the error is. That's good for outliers, but there's a problem...

def abs_gradient(diff):
    if diff > 0:
        return 1
    elif diff < 0:
        return -1
    else:
        return None  # undefined!

for d in [-2, -0.5, 0.5, 3]:
    analytical = abs_gradient(d)
    numerical  = numerical_gradient(abs_loss, d)
    print(f"diff={d:+.1f}: analytical={analytical:+}, numerical={numerical:+.6f}")

The gradient is undefined at diff = 0. The function $|x|$ has a sharp corner (a "kink") at $x = 0$: approaching from the left the slope is -1, approaching from the right the slope is +1. There is no single well-defined tangent line.

In mathematical terms: $|x|$ is not differentiable at $x = 0$.

At $|\text{diff}| = d$:

• Values: squared gives $d^2$, absolute gives $d$. These are only equal when $d = 1$ (a coincidence). For any other $d$, there's a jump in the function.
• Gradients: squared gives $2d$, absolute gives $1$. These are never equal (unless $d = 0.5$, another coincidence).

The function has discontinuities in both value and gradient at the boundary. Gradient descent can't work with jumps — we need smoothness.

With $d = 1$ the values happen to match (both give 1), but the slopes are visibly different — the squared curve arrives with a steep slope while the absolute line departs with slope 1. Try $d = 2$ to see a clear value jump too.

plt.figure(figsize=(8, 5))
plt.plot(diff, squared, '--', alpha=0.4, label='Squared (diff²)')
plt.plot(diff, absolute, '--', alpha=0.4, label='Absolute (|diff|)')
plt.plot(diff, naive, 'k-', linewidth=2, label='Naive combination')
plt.axvline(d, color='red', linestyle=':', alpha=0.5, label=f'±d = ±{d}')
plt.axvline(-d, color='red', linestyle=':', alpha=0.5)
plt.legend()
plt.xlabel('diff')
plt.ylabel('Loss')
plt.title('Naive combination — spot the jump at ±d')
plt.grid(True, alpha=0.3)
plt.show()

$A = 2d$, $B = -d^2$. The curve is now perfectly smooth at the transition point — no jump in value, no jump in gradient.

# From gradient matching: A = 2d
# From value matching:    d² = 2d·d + B  →  B = -d²

def smooth_combined_loss(diff, d):
    if abs(diff) <= d:
        return diff ** 2
    else:
        return 2 * d * abs(diff) - d ** 2

# Verify at boundary: both sides should give the same value and gradient
d = 1
print(f"Inner side at diff=d:  value = {d**2},  gradient = {2*d}")
print(f"Outer side at diff=d:  value = {2*d*d - d**2},  gradient = {2*d}")

# Plot
diff_range = np.linspace(-3, 3, 1000)
loss_values = [smooth_combined_loss(x, d) for x in diff_range]

plt.figure(figsize=(8, 5))
plt.plot(diff_range, loss_values, 'k-', linewidth=2, label='Smooth combined loss')
plt.plot(diff_range, diff_range ** 2, '--', alpha=0.3, label='diff²')
plt.plot(diff_range, 2*d*np.abs(diff_range) - d**2, '--', alpha=0.3, label=f'2d|diff| - d²')
plt.axvline(d, color='red', linestyle=':', alpha=0.5)
plt.axvline(-d, color='red', linestyle=':', alpha=0.5)
plt.legend()
plt.xlabel('diff')
plt.ylabel('Loss')
plt.title('Smooth combined loss — no jump, no kink!')
plt.grid(True, alpha=0.3)
plt.show()

def huber_loss(y_true, y_pred, delta):
    total = 0
    for y, yp in zip(y_true, y_pred):
        err = abs(y - yp)
        if err <= delta:
            total += 0.5 * (y - yp) ** 2
        else:
            total += delta * err - 0.5 * delta ** 2
    return total / len(y_true)

print(huber_loss([5, 2, 7], [4.8, 2.5, 10], delta=1))  # ≈ 0.8817
print(huber_loss([1, 2, 3], [1, 2, 3], delta=1))        # 0.0

The Huber curve hugs the MSE parabola near zero, then peels off and tracks the MAE line for large errors — smooth everywhere, no kinks, no explosions.

huber = np.where(
    np.abs(diff) <= delta,
    0.5 * diff ** 2,
    delta * np.abs(diff) - 0.5 * delta ** 2
)

plt.figure(figsize=(8, 5))
plt.plot(diff, mse, '--', alpha=0.6, label='½ diff² (MSE)')
plt.plot(diff, mae, '--', alpha=0.6, label='|diff| (MAE)')
plt.plot(diff, huber, 'k-', linewidth=2, label=f'Huber (δ={delta})')
plt.axvline(delta, color='red', linestyle=':', alpha=0.4)
plt.axvline(-delta, color='red', linestyle=':', alpha=0.4)
plt.legend(fontsize=11)
plt.xlabel('diff', fontsize=12)
plt.ylabel('Loss', fontsize=12)
plt.title('MSE vs MAE vs Huber Loss', fontsize=13)
plt.grid(True, alpha=0.3)
plt.ylim(-0.5, 6)
plt.show()

With the outlier, MSE explodes to ~1882 while Huber stays at ~19.8 — similar to MAE's 20.0 but with the advantage of being fully differentiable. On normal data, Huber behaves like MSE, giving smooth gradients near zero. Best of both worlds.

actual    = [1.0, 1.1, 0.9, 0.5, 3.0]
normal    = [0.9, 2.0, 0.7, 0.3, 2.0]
outlier   = [0.9, 2.0, 0.7, 0.3, 100.0]

for name, pred in [("Normal", normal), ("Outlier", outlier)]:
    print(f"{name} predictions:")
    print(f"  MSE:   {compute_mse(actual, pred):.4f}")
    print(f"  MAE:   {compute_mae(actual, pred):.4f}")
    print(f"  Huber: {huber_loss(actual, pred, delta=1):.4f}")
    print()