LAD Regression

lad_regression.py

import numpy as np
from cvxopt import matrix, solvers

def solve(X, y, l = 0.0):
  '''Solves the LAD regression problem, which can be formulated as:

      minimize || X a - y ||_1 + l * || a ||_1
          a

  X: numpy ndarray with shape (n, m)
  y: numpy array with shape (n,)
  l: regularization parameter. set to zero to perform unregularized optimization.

  returns the optimal solution a* as a numpy ndarray with shape (m,).'''

  n, m = X.shape
  assert y.shape == (n,)

  A = matrix(
      np.vstack([
        np.hstack([  X, -X, -np.eye(n) ]),
        np.hstack([ -X,  X, -np.eye(n) ]),
      ]))
  b = matrix(np.concatenate((y, -y)))
  c = matrix(np.concatenate((np.ones(2 * m) * l, np.ones(n))))
  solution = solvers.lp(c, A, b, solver = 'glpk')
  plus_and_minus = np.asarray(solution['x']).flatten()[:(2 * m)]
  return plus_and_minus[:m] - plus_and_minus[m:]

n = 1000
m = 10
X = np.random.randn(n, m) * (np.random.rand() + 3) + np.random.randn()
y = np.random.randn(n) * (np.random.rand() + 10) + np.random.randn()

a_star = solve(X, y)

print a_star
print np.sum(np.abs(np.dot(X, a_star) - y))
print np.sum(np.abs(np.dot(X, np.random.randn(*a_star.shape)) - y))
print np.sum(np.abs(np.dot(X, np.zeros_like(a_star)) - y))