Linear regression by Linear Programming

regression_lp.py

# (C) 2011 Mathieu Blondel
# License: BSD 3 clause

import numpy as np
import numpy.linalg as linalg
import pylab as pl
from cvxopt import matrix, solvers

np.random.seed(0)

def f(x):
    return 0.8 * x + 0.5

def add_noise(y, loc, scale):
    signs = np.sign(np.random.random(len(y)) - 0.5)
    return y + np.random.normal(loc, scale, size=len(y)) * signs

x = np.random.random(20)
y = add_noise(f(x), loc=0, scale=0.1)
indices = np.arange(len(x))
np.random.shuffle(indices)
indices = indices[:3]
y[indices] = add_noise(y[indices], loc=0, scale=1.0)

# least-squares solution
# min ||X^T w - y||^2 where w = [beta, alpha]

X = np.array([x, np.ones(len(x))]).T
w_ls = linalg.solve(np.dot(X.T, X), np.dot(X.T, y))

def f_ls(x):
    return x * w_ls[0] + w_ls[1]

# linear programming solution
# min \sum_i e_i
# s.t. -e_i <= alpha + beta x_i - y_i <= e_i

c = np.ones(len(x) + 2)
c[-2:] = 0

G11 = G12 = np.diag(-np.ones(len(x)))
G21 = np.array([-x, -np.ones(len(x))])
G22 = np.array([+x, np.ones(len(x))])

G_left = np.vstack((G11, G21))
G_right = np.vstack((G12, G22))
G = np.hstack((G_left, G_right))

h = np.ravel(np.vstack((-y, y)))

c = matrix(c.tolist())
G = matrix(G.tolist())
h = matrix(h.tolist())
sol = solvers.lp(c, G, h)
w_lp = list(sol['x'])[-2:]

def f_lp(x):
    return x * w_lp[0] + w_lp[1]

pl.plot([-0.5, 1.5], [f_ls(-0.5), f_ls(1.5)], "r--", label="least-squares")
pl.plot([-0.5, 1.5], [f_lp(-0.5), f_lp(1.5)], "b--", label="linear programming")
pl.scatter(x, y, color="black", label="training points")

pl.legend(loc='upper left')

pl.show()