Huber Loss & Smooth L1

1
#include <iostream>
2
#include <cmath>
3
#include <vector>
4
5
// Huber loss with parameter delta
6
double huber_loss(double r, double delta) {
7
    double ar = std::fabs(r);
8
    if (ar <= delta) return 0.5 * r * r; // quadratic region
9
    return delta * ar - 0.5 * delta * delta; // linear region with continuity
10
}
11
12
// Derivative of Huber w.r.t. residual r
13
// Note: gradient w.r.t. prediction y_hat is -huber_grad(r, delta)
14
double huber_grad(double r, double delta) {
15
    double ar = std::fabs(r);
16
    if (ar <= delta) return r; // slope grows with r
17
    return delta * (r < 0 ? -1.0 : 1.0); // saturated slope
18
}
19
20
// Smooth L1 (scaled Huber) with parameter beta
21
// s_beta(r) = (1/beta) * huber_beta(r)
22
double smooth_l1(double r, double beta) {
23
    double ar = std::fabs(r);
24
    if (ar < beta) return 0.5 * r * r / beta; // quadratic with 1/beta factor
25
    return ar - 0.5 * beta; // linear tail
26
}
27
28
// Derivative of Smooth L1 w.r.t. residual r
29
double smooth_l1_grad(double r, double beta) {
30
    double ar = std::fabs(r);
31
    if (ar < beta) return r / beta; // continuous slope 1 at boundary
32
    return (r < 0 ? -1.0 : 1.0);
33
}
34
35
int main() {
36
    std::vector<double> residuals = { -3.0, -0.2, 0.0, 0.1, 2.5 };
37
    double delta = 1.0; // Huber threshold
38
    double beta  = 1.0; // Smooth L1 threshold
39
40
    std::cout << "r\tHuber\tHuber'\tSmoothL1\tSmoothL1'\n";
41
    for (double r : residuals) {
42
        double hl = huber_loss(r, delta);
43
        double hg = huber_grad(r, delta);
44
        double sl = smooth_l1(r, beta);
45
        double sg = smooth_l1_grad(r, beta);
46
        std::cout << r << '\t' << hl << '\t' << hg << '\t' << sl << '\t' << sg << '\n';
47
    }
48
    return 0;
49
}
50

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
struct DataPoint { double x, y; };
5
6
// Huber gradient w.r.t residual r
7
double huber_grad(double r, double delta) {
8
    double ar = fabs(r);
9
    if (ar <= delta) return r;
10
    return delta * (r < 0 ? -1.0 : 1.0);
11
}
12
13
int main() {
14
    // Synthetic data: y = 2x + 1 with one outlier
15
    vector<DataPoint> data;
16
    for (int i = 0; i <= 10; ++i) {
17
        double x = i / 2.0;
18
        double y = 2.0 * x + 1.0 + 0.05 * ((i % 3) - 1); // small noise
19
        data.push_back({x, y});
20
    }
21
    // Inject a strong outlier
22
    data.push_back({3.5, 30.0});
23
24
    // Initialize parameters for line y = a x + b
25
    double a_mse = 0.0, b_mse = 0.0;
26
    double a_hub = 0.0, b_hub = 0.0;
27
28
    double lr = 0.05;     // learning rate
29
    double delta = 1.0;   // Huber threshold
30
    int iters = 400;
31
32
    auto step_mse = [&](double &a, double &b) {
33
        double ga = 0.0, gb = 0.0; // gradients
34
        for (auto &p : data) {
35
            double yhat = a * p.x + b;
36
            double r = p.y - yhat; // residual
37
            // MSE loss: 0.5 r^2 -> dL/dyhat = -r
38
            ga += -r * p.x;
39
            gb += -r;
40
        }
41
        ga /= data.size(); gb /= data.size();
42
        a -= lr * ga; b -= lr * gb;
43
    };
44
45
    auto step_huber = [&](double &a, double &b) {
46
        double ga = 0.0, gb = 0.0; // gradients
47
        for (auto &p : data) {
48
            double yhat = a * p.x + b;
49
            double r = p.y - yhat; // residual
50
            double g = huber_grad(r, delta); // dL/dr
51
            // chain rule: dL/dyhat = -g
52
            ga += -g * p.x;
53
            gb += -g;
54
        }
55
        ga /= data.size(); gb /= data.size();
56
        a -= lr * ga; b -= lr * gb;
57
    };
58
59
    for (int t = 0; t < iters; ++t) {
60
        step_mse(a_mse, b_mse);
61
        step_huber(a_hub, b_hub);
62
    }
63
64
    cout.setf(std::ios::fixed); cout << setprecision(4);
65
    cout << "Ground truth: a=2.0000 b=1.0000\n";
66
    cout << "With outlier -> MSE fit:   a=" << a_mse << " b=" << b_mse << "\n";
67
    cout << "With outlier -> Huber fit: a=" << a_hub << " b=" << b_hub << "\n";
68
69
    return 0;
70
}
71

1
#include <iostream>
2
#include <vector>
3
#include <cmath>
4
#include <iomanip>
5
6
// Compute Smooth L1 loss and gradient per coordinate for vectors
7
// y_pred and y_true must be same length; beta is per-dimension or scalar
8
9
struct SmoothL1Result {
10
    double loss;                 // summed loss
11
    std::vector<double> grad;    // gradient w.r.t. y_pred
12
};
13
14
SmoothL1Result smooth_l1_vec(const std::vector<double>& y_pred,
15
                             const std::vector<double>& y_true,
16
                             const std::vector<double>& beta) {
17
    size_t d = y_pred.size();
18
    SmoothL1Result res; res.loss = 0.0; res.grad.assign(d, 0.0);
19
    for (size_t j = 0; j < d; ++j) {
20
        double r = y_true[j] - y_pred[j];
21
        double b = beta[j];
22
        double ar = std::fabs(r);
23
        if (ar < b) {
24
            res.loss += 0.5 * r * r / b;
25
            res.grad[j] = -(r / b); // dL/dyhat = - dL/dr
26
        } else {
27
            res.loss += ar - 0.5 * b;
28
            res.grad[j] = -(r < 0 ? -1.0 : 1.0);
29
        }
30
    }
31
    return res;
32
}
33
34
int main() {
35
    // Example: bounding boxes as (cx, cy, w, h)
36
    std::vector<double> y_true = {50.0, 40.0, 120.0, 80.0};
37
    std::vector<double> y_pred = {52.5, 38.0, 130.0, 70.0};
38
39
    // Per-dimension beta (common to set smaller beta for widths/heights)
40
    std::vector<double> beta   = {1.0, 1.0, 1.0, 1.0};
41
42
    SmoothL1Result r = smooth_l1_vec(y_pred, y_true, beta);
43
44
    std::cout << std::fixed << std::setprecision(4);
45
    std::cout << "Smooth L1 loss (sum over dims): " << r.loss << "\nGradients w.r.t. y_pred:" << std::endl;
46
    for (double g : r.grad) std::cout << g << ' ';
47
    std::cout << std::endl;
48
49
    return 0;
50
}
51