L1标准而不是L2规范

Question

i我想知道Python中是否有一个与

scipy.linalg.lstsq

的功能相同的工作，但使用“最小绝对偏差”回归而不是“最小二乘”回归（OLS）。我想使用

L1

规范，而不是

L2

规范。

实际上，我有3D点，我想要它们最合适的平面。常见的方法是从最少的正方形方法，例如这种github

link。但是，众所周知，这并不总是能带来最佳拟合度，尤其是当我们在数据集中有闯入者时。最好计算最小绝对偏差。两种方法之间的差异被更多地解释了。 由于它是一个Ax = b矩阵方程，因此不会通过诸如MAD之类的函数来求解，并且需要循环以最大程度地减少结果。我想知道是否有人知道Python中的相关功能 - 可能是在线性代数软件包中 - 可以计算“最小绝对偏差”回归？

使用

scipy.optimize.minimize

和自定义

cost_function

Answer 1

我们首先导入必需品，

from scipy.optimize import minimize
import numpy as np

定义自定义成本函数（以及用于获取拟合值的便利包装），

def fit(X, params):
    return X.dot(params)


def cost_function(params, X, y):
    return np.sum(np.abs(y - fit(X, params)))

，然后，如果您有一些

X

（设计矩阵）和

y

（观察），我们可以做以下操作，

output = minimize(cost_function, x0, args=(X, y))

y_hat = fit(X, output.x)

where

x0

是对最佳参数的一些合适的初始猜测（您可以在此处获取@Jamesphillips的建议，并使用OLS方法中的拟合参数）。

在任何情况下，当测试以某种人为的例子进行测试时，

X = np.asarray([np.ones((100,)), np.arange(0, 100)]).T
y = 10 + 5 * np.arange(0, 100) + 25 * np.random.random((100,))

我发现，

      fun: 629.4950595335436
 hess_inv: array([[  9.35213468e-03,  -1.66803210e-04],
       [ -1.66803210e-04,   1.24831279e-05]])
      jac: array([  0.00000000e+00,  -1.52587891e-05])
  message: 'Optimization terminated successfully.'
     nfev: 144
      nit: 11
     njev: 36
   status: 0
  success: True
        x: array([ 19.71326758,   5.07035192])

和，

fig = plt.figure()
ax = plt.axes()

ax.plot(y, 'o', color='black')
ax.plot(y_hat, 'o', color='blue')

plt.show()

以蓝色的拟合值和黑色的数据。

您可以使用scipy.minimimize函数解决问题。您必须设置要最小化的功能（在我们的情况下，具有Z = AX + By + C的平面）和错误函数（L1 Norm），然后以一些起始值运行最小化器。

import numpy as np
import scipy.linalg
from scipy.optimize import minimize
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

def fit(X, params):
    # 3d Plane Z = aX + bY + c
    return X.dot(params[:2]) + params[2]

def cost_function(params, X, y):
    # L1- norm
    return np.sum(np.abs(y - fit(X, params)))

我们生成3D点

# Generating  3-dim points
mean = np.array([0.0,0.0,0.0])
cov = np.array([[1.0,-0.5,0.8], [-0.5,1.1,0.0], [0.8,0.0,1.0]])
data = np.random.multivariate_normal(mean, cov, 50)

Answer 2

output = minimize(cost_function, [0.5,0.5,0.5], args=(np.c_[data[:,0], data[:,1]], data[:, 2]))
y_hat = fit(np.c_[data[:,0], data[:,1]], output.x)

X,Y = np.meshgrid(np.arange(min(data[:,0]), max(data[:,0]), 0.5),    np.arange(min(data[:,1]), max(data[:,1]), 0.5))
XX = X.flatten()
YY = Y.flatten()


# # evaluate it on grid
Z = output.x[0]*X + output.x[1]*Y + output.x[2]
fig = plt.figure(figsize=(10,10))
ax = fig.gca(projection='3d')
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, alpha=0.2)
ax.scatter(data[:,0], data[:,1], data[:,2], c='r')
plt.show()

注：我已经使用了先前的响应代码，而从github使用了

代码

基于公式

https：//home.cs.colorado.edu/~srirams/courses/csci5654-fall13/regress-sli.pdf 我创建了这个Python脚本。 对于较大的数据集，不幸的是，它变得非常慢。

import numpy as np
import scipy as sp
from scipy.optimize import linprog
from scipy.sparse import csc_matrix, hstack




# Create matrix A and vector y
A = motor_filtered[['x1','x2']].values
y = motor_filtered['y'].values

#l1 problem using https://home.cs.colorado.edu/~srirams/courses/csci5654-fall13/regression-sli.pdf

A_sparse = csc_matrix(A)

x_dim, y_dim = A.shape[1], y.shape[0]
eye = np.eye(x_dim)

zeros = np.zeros((1, x_dim))
ones = np.ones((1, A_sparse.shape[0]))

# Concatenate the sparse vectors
obj = np.hstack((zeros, ones))


I = sp.sparse.eye(A_sparse.shape[0])
lhs_ineq = sp.sparse.hstack((A_sparse, -I))


lhs_ineq = sp.sparse.vstack((lhs_ineq, sp.sparse.hstack((-A_sparse, -I))))


rhs_ineq =  np.hstack((y,-y))

bnd = [*((None, None) for _ in range(x_dim)), *((0, None) for _ in range(A_sparse.shape[0]))]

res = linprog(c=obj, A_ub=lhs_ineq, b_ub=rhs_ineq,  bounds=bnd, method="highs")
print(res)
print(res.x[:x_dim])