Huber Fitting

Consider the robust regression problem

\[ \text{minimize} ~\sum_{i=1}^m\text{huber}(a_i^T x - b_i), \]

where

\[ \text{huber}(x) = \left\{\begin{aligned} &(1/2)x^2 & |x| \leq 1 \\\ &|x| - (1/2) & |x| > 1 \end{aligned} \right. \]

which has the graph form representation

\[ \begin{aligned} &\text{minimize} & & \sum_{i=1}^m\text{huber}(y_i - b_i) \\
& \text{subject to} & & y = A x. \end{aligned} \]

or equivalently

\[ \begin{aligned} &\text{minimize} & & f(y) + g(x) \\
& \text{subject to} & & y = A x, \end{aligned} \]

where

\[ f_i(y_i) = \text{huber}(y_i - b_i), ~~\text{ and } ~~g_j(x_j) = 0. \]

MATLAB Code

1
2
3
4
5
6
7
8
9
10
11
% Generate Data
A = randn(100, 10);
b = randn(100, 1);

% Populate f and g
f.h = kHuber;
f.b = b;
g.h = kZero;

% Solve
x = pogs(A, f, g);

A similar example can be found in the file 

1
<pogs>/examples/matlab/huber_fit.m
.

R Code

1
2
3
4
5
6
7
8
9
10
# Generate Data
A = matrix(rnorm(100 * 10), 100, 10)
b = rnorm(100)

# Populate f and g
f = list(h = kHuber(), b = b)
g = list(h = kZero())

# Solve
solution = pogs(A, f, g)

A similar example can be found in the file 

1
<pogs>/examples/r/huber_fit.R
.

C++ Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
#include <random>
#include <vector>

#include "pogs.h"

int main() {
  // Generate Data
  size_t m = 100, n = 10;
  std::vector<double> A(m * n);
  std::vector<double> b(m);
  std::vector<double> x(n);
  std::vector<double> y(m);

  std::default_random_engine generator;
  std::normal_distribution<double> n_dist(0.0, 1.0);

  for (unsigned int i = 0; i < m * n; ++i)
    A[i] = n_dist(generator);

  for (unsigned int j = 0; j < n; ++j)
    b[i] = n_dist(generator);

  // Populate f and g
  PogsData<double, double*> pogs_data(A.data(), m, n);
  pogs_data.x = x.data();
  pogs_data.y = y.data();

  pogs_data.f.reserve(m);
  for (unsigned int i = 0; i < m; ++i)
    pogs_data.f.emplace_back(kHuber, 1.0, b[i]);

  pogs_data.g.reserve(n);
  for (unsigned int i = 0; i < n; ++i)
    pogs_data.g.emplace_back(kZero);

  // Solve
  Pogs(&pogs_data);
}

A similar example can be found in the file 

1
<pogs>/examples/cpp/huber_fit.cpp
.