Huber regression¶

$$ \newcommand{\eg}{{\it e.g.}} \newcommand{\ie}{{\it i.e.}} \newcommand{\argmin}{\operatornamewithlimits{argmin}} \newcommand{\mc}{\mathcal} \newcommand{\mb}{\mathbb} \newcommand{\mf}{\mathbf} \newcommand{\minimize}{{\text{minimize}}} \newcommand{\argmin}{{\text{argmin}}} \newcommand{\diag}{{\text{diag}}} \newcommand{\cond}{{\text{cond}}} \newcommand{\rank}{{\text{rank }}} \newcommand{\range}{{\mathcal{R}}} \newcommand{\null}{{\mathcal{N}}} \newcommand{\tr}{{\text{trace}}} \newcommand{\dom}{{\text{dom}}} \newcommand{\dist}{{\text{dist}}} \newcommand{\E}{\mathbf{E}} \newcommand{\var}{\mathbf{var}} \newcommand{\R}{\mathbf{R}} \newcommand{\B}{\mathbf{B}} \newcommand{\SM}{\mathbf{S}} \newcommand{\ball}{\mathcal{B}} \newcommand{\onev}{\mathbf{1}} \newcommand{\bmat}[1]{\begin{bmatrix}#1\end{bmatrix}} $$

Jong-Han Kim (jonghank@inha.ac.kr)

---

Standard regression¶

In a standard regression problem we are given data $(x_i,y_i)\in \R^d \times \R$, $i=1,\ldots, n$, and fit a linear (affine) model

$$\hat y_i = \theta ^T x_i + b$$

where $\theta \in\R^d$ and $b \in {\bf R}$.

The residuals are $r_i = \hat y_i - y_i$. In standard (least-squares) regression we choose $\theta$ and $b$ to minimize $\|r\|_2^2 = \sum_i r_i^2$. For this choice of $\theta$ and $b$ the mean of the optimal residuals is zero.

In a vector form, we can write the problem as

$$ \hat{y} = \bmat{X & \onev}\bmat{\theta \\ b} $$

where

$$ X = \bmat{x_1^T \\ x_2^T \\ \vdots \\ x_n^T}\in\R^{n\times d} $$

So the problem is to find $\theta\in\R^d$ and $b\in\R$ for

$$ \begin{aligned} \underset{\theta, b}{\minimize} \quad & \|X\theta + \onev b - y\|^2 \end{aligned} $$

A simple variant is to add (Tikhonov) regularization, meaning we solve the optimization problem

$$\begin{array}{ll} \mbox{minimize} & \|X\theta + \onev b - y\|_2^2 + \lambda \|\theta \|_2^2, \end{array}$$

with some $\lambda>0$.

---

Robust (Huber) regression¶

A more sophisticated variant is to replace the square function, $\|X\theta + \onev b - y\|_2^2$, with the Huber function

$$ h_\rho(x) = \begin{cases} x^2 &\quad \text{if } |x| \le \rho \\ \rho\left(2|x|-\rho\right) &\quad \text{if }|x| \gt \rho \end{cases} $$

where $\rho>0$ is the Huber threshold. The image below shows the square function on the left and the Huber function on the right.

Huber regression is the same as standard (least-squares) regression for small residuals, but allows (some) large residuals.

---

Numerical example¶

In the following code we do a numerical example of Huber regression.

We generate $n=450$ measurements with $d=300$ feasures.

We randomly choose $\theta^\text{true}$, $b^\text{true}$, and $x_i \sim \mc N(0,I)$.

We set $y_i = (\theta^\text{true})^Tx_i + b^\text{true} + \epsilon_i$, where $\epsilon_i \sim \mc N(0,1)$ represents the measurement noise.

Then all kinds of nasty things come in; we replace $y_i$ with $-y_i$ with probability $p$.

The data has fraction $p$ of (non-obvious) wrong measurements.The distribution of "good" and "bad" $y_i$ are the same.

Our goal is to recover $\theta^\text{true} \in\R^d$ and $b^\text{true} \in\R$ from the measurements $y\in \R^n$. We compare three approaches:

• Standard regression (least squares)
• Huber regression
• "Prescient" regression, where we know which measurements had their sign flipped.

Note that the "prescient" approach should be equivalent to "no-flip" case ($p=0$), since knowing the flipped measurements implies that we can re-flip those to recover the original unflipped signals. So the perfornance of the "prescient" approaches with different $p$'s should be identical for all $p$'s.

We generate $50$ problem instances, with $p$ varying from $0$ to $0.15$, and plot the relative error in reconstructing $\theta^\text{true}$ for the three approaches.

Notice that in the range $p \in [0,0.06]$, Huber regression almost matches prescient regression. Standard regression, by contrast, fails even for very small $p$.