What: I'm trying to minimise the Huber-Loss for a linear regression using Stochastic Gradient Descent from scratch.
Problem: It seems like that the coeffcient $m$ doesn't get optimised, therefore the fitted line is nowhere near the OLS regression line which serves as reference.
Expectation: Optimised coefficient resulting in a fitted regression line close to the reference OLS regression line.
Explanation: The Huber-Loss is defined as
$$ L_{\delta}(a) = \begin{array}{l}\frac{1}{2}a^2 \quad \textrm{for} \quad \vert a \vert \leq \delta \\
\delta (\vert a \vert - \frac{1}{2}\delta) \quad \textrm{otherwise} \end{array} $$
The variable $a$ refers to the residual $a = y - \hat{y}$, which means for linear regression: $a = y - (mx+n)$. To use Stochastic Gradient Descent, the partial derivatives with respect to the coefficient $m$ and the intercept $n$ have to be calculated:
$$\nabla L_{\delta}(m, n, y) = \left[\begin{array}{l}
\dfrac{\partial L}{\partial m} = \begin{array}{l} -x\ (-mx+n-y)\quad \textrm{for} \quad \vert y-(mx+n) \vert \leq \delta \\
\frac{\delta x(mx+n-y)}{\vert mx+n-y \vert} \quad \textrm{otherwise}
\end{array}\\
\dfrac{\partial L}{\partial n} = \begin{array}{l} mx+n-y\quad \textrm{for} \quad \vert y-(mx+n) \vert \leq \delta \\
\frac{\delta (mx+n-y)}{\vert mx+n-y \vert} \quad \textrm{otherwise}
\end{array} \\
\end{array}\right]$$
I implemented the above in R the following way.
grad_m_huber <- function(
Y,
X,
m,
n,
delta
) {
loss <- ifelse(
abs(Y - (m * X + n)) <= delta,
-X * (-m * X - n + Y),
(delta * X * (m * X + n - Y)) / abs(n + m * X - Y)
)
return(loss)
}
grad_n_huber <- function(
Y,
X,
m,
n,
delta
) {
loss <- ifelse(
abs((m * X + n) - Y) <= delta,
m * X + n - Y,
(delta * (m * X + n - Y)) / abs(m * X + n - Y)
)
return(loss)
}
huber_sgd <- function(
X,
Y,
epochs,
lr,
batch_size,
delta
) {
m <- 0
n <- 0
for(i in 1:epochs) {
batch <- sample(length(X), batch_size)
Y_batch <- Y[batch]
X_batch <- X[batch]
g_m <- sum(grad_m_huber(Y_batch, X_batch, m, n, delta))
g_n <- sum(grad_n_huber(Y_batch, X_batch, m, n, delta))
m <- m - lr * g_m
n <- n - lr * g_n
}
return(
list(
grad_m = g_m,
grad_n = g_n,
m = m,
n = n
)
)
}
But when I run the Huber-Loss regression with a delta close to 0, which should result in a fitting line close to a linear regression (ols_mod <- lm(Y ~ X)):
set.seed(303)
Y <- rnorm(550, 155, 15)
scale_factor <- .27
eps <- rnorm(550, 15, 5)
X <- (Y * scale_factor) + eps
huber_sgd_mod <- huber_sgd(
X,
Y,
epochs = 5000,
epsillon = .001,
batch_size = 1,
delta = .01
)
I get
As the plotted loss indicated that SGD didn't converge yet I changed the number of epochs to epochs = 10000 and ran it again. It seems like that only the intercept gets optimised.
I also compared my from scratch solution to a Huber-Regression calculated by MASS::rlm(Y ~ X, psi = MASS::psi.huber) which shows the expected fitted line close to the OLS regression line.
Does anyone know what am I doing wrong here?
