Huber fitting

Huber fitting or the robust least-squares problem performs linear regression under the assumption that there are outliers in the data. The fitting problem is written as

\[\begin{array}{ll} \mbox{minimize} & \sum_{i=1}^{m} \phi_{\rm hub}(a_i^T x - b_i), \end{array}\]

with the Huber penalty function \(\phi_{\rm hub}:\mathbf{R}\to\mathbf{R}\) defined as

\[\begin{split}\phi_{\rm hub}(u) = \begin{cases} u^2 & |u| \le 1 \\ (2|u| - 1) & |u| > 1 \end{cases}\end{split}\]

The problem has the following equivalent form,

\[\begin{split}\begin{array}{ll} \mbox{minimize} & \frac{1}{2} u^T u + \boldsymbol{1}^T v \\ \mbox{subject to} & -u-v \le Ax - b \le u+v \\ & 0 \le u \le 1 \\ & v \ge 0 \end{array}\end{split}\]

Python

import osqp
import numpy as np
import scipy as sp
import scipy.sparse as sparse

# Generate problem data
sp.random.seed(1)
n = 10
m = 100
Ad = sparse.random(m, n, density=0.5, format='csc')
x_true = np.random.randn(n) / np.sqrt(n)
ind95 = (np.random.rand(m) < 0.95).astype(float)
b = Ad.dot(x_true) + np.multiply(0.5*np.random.randn(m), ind95) \
    + np.multiply(10.*np.random.rand(m), 1. - ind95)

# OSQP data
Im = sparse.eye(m)
P = sparse.block_diag((sparse.csc_matrix((n, n)), Im,
                       sparse.csc_matrix((m, m))), format='csc')
q = np.append(np.zeros(m+n), np.ones(m))
A = sparse.vstack([
        sparse.hstack([Ad, Im, Im]),
        sparse.hstack([Ad, -Im, -Im]),
        sparse.hstack([sparse.csc_matrix((2*m, n)), sparse.eye(2*m)])
        ]).tocsc()
l = np.hstack([b, -np.inf*np.ones(m), np.zeros(2*m)])
u = np.hstack([np.inf*np.ones(m), b, np.ones(m), np.inf*np.ones(m)])

# Create an OSQP object
prob = osqp.OSQP()

# Setup workspace
prob.setup(P, q, A, l, u)

# Solve problem
res = prob.solve()

Matlab

% Generate problem data
rng(1)
n = 10;
m = 100;
Ad = sprandn(m, n, 0.5);
x_true = randn(n, 1) / sqrt(n);
ind95 = rand(m, 1) > 0.95;
b = Ad*x_true + 10*rand(m, 1).*ind95 + 0.5*randn(m, 1).*(1-ind95);

% OSQP data
Im = speye(m);
P = blkdiag(sparse(n, n), Im, sparse(m, m));
q = [zeros(m + n, 1); ones(m, 1)];
A = [Ad, Im, Im;
     Ad, -Im, -Im;
     sparse(2*m, n), speye(2*m)];
l = [b; -inf*ones(m, 1); zeros(2*m, 1)];
u = [inf*ones(m, 1); b; ones(m, 1); inf*ones(m, 1)];

% Create an OSQP object
prob = osqp;

% Setup workspace
prob.setup(P, q, A, l, u);

% Solve problem
res = prob.solve();