# Efficient method for Laplace regression

I want to calculate numerically the maximum likelihood estimators of $(\beta,\sigma)$ for the linear regression model:

$$y_j = x_j^{\top}\beta + \epsilon_j,$$

where $j=1,\dots,n$, $\beta$ is $p$-dimensional, and $\epsilon_j$ are i.i.d. according to a Laplace distribution with location zero and scale $\sigma$.

Given the differentiability issues with this distribution, I cannot obtain the MLE of $\beta$ by solving the score functions. I need to do this for a number of data sets with large $p$. Is there an efficient method for doing so?

The MLE minimizes

$\min \| X \beta - y\|_{1}$

Unfortunately, there isn't any simple closed form solution to this optimization problem. However, this is a convex optimization problem for which there are many available approaches.

For reasonably small instances (e.g. $n$ on the order of $100,000$ and $p$ is less than say $1,000$), the simplest approach is to use a linear programming library routine to solve the linear programming problem:

$\min \sum_{i=1}^{n} t_{i}$

subject to

$t \geq X\beta-y$

$t \geq y-X\beta$.

For larger instances (e.g. $n$ is in the billions or even larger), you might consider a stochastic subgradient descent method. This is relatively easy to implement in a "big data" environment with hadoop.

• Many thanks for your answer. Your suggestions look very interesting. – TheRock Apr 8 '15 at 14:26