约束线性回归/二次规划python答案

【问题标题】：constrained linear regression / quadratic programming python约束线性回归/二次规划python
【发布时间】：2016-10-04 12:55:29
【问题描述】：

我有一个这样的数据集：

import numpy as np

a = np.array([1.2, 2.3, 4.2])
b = np.array([1, 5, 6])
c = np.array([5.4, 6.2, 1.9])

m = np.vstack([a,b,c])
y = np.array([5.3, 0.9, 5.6])

并且想要拟合一个受约束的线性回归

y = b1*a + b2*b + b3*c

其中所有 b 的和为 1 且为正数：b1+b2+b3=1

此处指定了 R 中的一个类似问题：

https://stats.stackexchange.com/questions/21565/how-do-i-fit-a-constrained-regression-in-r-so-that-coefficients-total-1

如何在 python 中做到这一点？

【问题讨论】：

标签： python scipy linear-regression quadratic-programming

【解决方案1】：

编辑： 这两种方法非常通用，适用于中小型实例。要获得更有效的方法，请检查 chthonicdaemon 的答案（使用自定义预处理和 scipy 的 optimize.nnls）。

使用 scipy

代码

import numpy as np
from scipy.optimize import minimize

a = np.array([1.2, 2.3, 4.2])
b = np.array([1, 5, 6])
c = np.array([5.4, 6.2, 1.9])

m = np.vstack([a,b,c])
y = np.array([5.3, 0.9, 5.6])

def loss(x):
    return np.sum(np.square((np.dot(x, m) - y)))

cons = ({'type': 'eq',
         'fun' : lambda x: np.sum(x) - 1.0})

x0 = np.zeros(m.shape[0])
res = minimize(loss, x0, method='SLSQP', constraints=cons,
               bounds=[(0, np.inf) for i in range(m.shape[0])], options={'disp': True})

print(res.x)
print(np.dot(res.x, m.T))
print(np.sum(np.square(np.dot(res.x, m) - y)))

输出

Optimization terminated successfully.    (Exit mode 0)
        Current function value: 18.817792344
        Iterations: 5
        Function evaluations: 26
        Gradient evaluations: 5
[ 0.7760881  0.         0.2239119]
[ 1.87173571  2.11955951  4.61630834]
18.817792344

评估

很明显，模型能力/模型复杂度不足以获得良好的性能（高损失！）

使用由 cvxpy 建模的通用 QP/SOCP 优化

优点：

cvxpy 证明问题是凸的
保证了凸优化问题的收敛性（上述情况可能也是如此）
总的来说：准确度更高
总的来说：在数值不稳定性方面更稳健（求解器只能求解 SOCP；而不是像上面的 SLSQP 方法这样的非凸模型！）

代码

import numpy as np
from cvxpy import *

a = np.array([1.2, 2.3, 4.2])
b = np.array([1, 5, 6])
c = np.array([5.4, 6.2, 1.9])

m = np.vstack([a,b,c])
y = np.array([5.3, 0.9, 5.6])

X = Variable(m.shape[0])
constraints = [X >= 0, sum_entries(X) == 1.0]

product = m.T * diag(X)
diff = sum_entries(product, axis=1) - y
problem = Problem(Minimize(norm(diff)), constraints)
problem.solve(verbose=True)

print(problem.value)
print(X.value)

输出

ECOS 2.0.4 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS

It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT
 0  +0.000e+00  -0.000e+00  +2e+01  5e-01  1e-01  1e+00  4e+00    ---    ---    1  1  - |  -  - 
 1  +2.451e+00  +2.539e+00  +4e+00  1e-01  2e-02  2e-01  8e-01  0.8419  4e-02   2  2  2 |  0  0
 2  +4.301e+00  +4.306e+00  +2e-01  5e-03  7e-04  1e-02  4e-02  0.9619  1e-02   2  2  2 |  0  0
 3  +4.333e+00  +4.334e+00  +2e-02  4e-04  6e-05  1e-03  4e-03  0.9326  2e-02   2  1  2 |  0  0
 4  +4.338e+00  +4.338e+00  +5e-04  1e-05  2e-06  4e-05  1e-04  0.9698  1e-04   2  1  1 |  0  0
 5  +4.338e+00  +4.338e+00  +3e-05  8e-07  1e-07  3e-06  7e-06  0.9402  7e-03   2  1  1 |  0  0
 6  +4.338e+00  +4.338e+00  +7e-07  2e-08  2e-09  6e-08  2e-07  0.9796  1e-03   2  1  1 |  0  0
 7  +4.338e+00  +4.338e+00  +1e-07  3e-09  4e-10  1e-08  3e-08  0.8458  2e-02   2  1  1 |  0  0
 8  +4.338e+00  +4.338e+00  +7e-09  2e-10  2e-11  9e-10  2e-09  0.9839  5e-02   1  1  1 |  0  0

OPTIMAL (within feastol=1.7e-10, reltol=1.5e-09, abstol=6.5e-09).
Runtime: 0.000555 seconds.

4.337947939  # needs to be squared to be compared to scipy's output!
             #  as we are using l2-norm (outer sqrt) instead of sum-of-squares
             #  which is nicely converted to SOCP-form and easier to
             #  tackle by SOCP-based solvers like ECOS
             #  -> does not change the solution-vector x, only the obj-value
[[  7.76094262e-01]
 [  7.39698388e-10]
 [  2.23905737e-01]]

【讨论】：

这很棒。太感谢了。我只是好奇为什么在最小化函数中调用损失函数不需要括号中的参数？我们将损失函数定义为 loss(x)，但是当我们在最小化函数中调用它时，我们只输入 loss，而不是 loss(x)。你知道为什么吗？？
@sascha 这是一个有用的解释！请注意，在cvxpy 功能中，我必须更改sum_entries 对sum 的调用——我这样做没有 import *，因为它的价值。

【解决方案2】：

你可以通过一点数学和scipy.optimize.nnls得到一个很好的解决方案：

首先我们做数学：

如果

y = b1*a + b2*b + b3*c 和 b1 + b2 + b3 = 1，然后 b3 = 1 - b1 - b2。

如果我们替换并简化，我们最终会得到 p>

y - c = b1(a - c) + b2(b - c)

现在，我们没有任何等式约束，nnls 可以直接求解：

import scipy.optimize
A = np.vstack([a - c, b - c]).T
(b1, b2), norm = scipy.optimize.nnls(A, y - c)
b3 = 1 - b1 - b2

这恢复了与使用 cvxpy 在另一个答案中获得的相同的解决方案。

b1 = 0.77608809648662802
b2 = 0.0
b3 = 0.22391190351337198
norm = 4.337947941595865

这种方法可以推广到任意数量的维度，如下所示。假设我们有一个矩阵 B，矩阵 B 由原始问题中的 a、b、c 构成，列中排列。任何其他维度都将添加到此。

现在，我们可以做

A = B[:, :-1] - B[:, -1:]
bb, norm = scipy.optimize.nnls(A, y - B[:, -1])
bi = np.append(bb, 1 - sum(bb))

【讨论】：

这是一个不错的方法！
@sascha 谢谢！我很想看看 nnls 与更通用的求解器相比如何解决更大的问题。
对于固定数量的行，this fairly naive benchmark 表明您的方法非常快。但是如果让它更通用地解决 x 行的问题会很有趣。使用更多行，我的基于 SOCP 的方法也比 SLSQP 更快，后者占 3 行。

【解决方案3】：

关于sascha's scipy implementation 的一条评论：请注意，使用 scipy 最小化，SLSQP 的反复试验性质可能会为您提供一个略微的解决方案，除非您制定一些其他规范，即最大迭代次数 (maxiter) 和最大容差 (ftol)，详见 scipy 文档here。

默认值为：maxiter=100 和 ftol=1e-06。

这里有一个例子来说明使用矩阵表示法：首先摆脱约束和界限。为简单起见，还假设截距 = 0。在这种情况下，任何多元回归的系数，如第 4 页 here 所述，将（准确地说）：

def betas(y, x):
    # y and x are ndarrays--response & design matrixes
    return np.dot(np.linalg.inv(np.dot(x.T, x)), np.dot(x.T, y))

现在，鉴于最小二乘回归的目标是最小化残差平方和，采用 sascha 的损失函数（稍微重写）：

def resids(b, y, x):
    resid = y - np.dot(x, b)
    return np.dot(resid.T, resid)

根据您的实际 Y 和 X 向量，您可以将上述第一个等式中生成的“真实”β 代入第二个等式，以获得更好的“基准”。将此基准与 res 的 .fun 属性（scipy最小化吐出的内容）进行比较。即使是微小的变化也会导致结果系数发生有意义的变化。

所以长话短说，使用类似的东西会牺牲速度但会提高准确性

options={'maxiter' : 1000, 'ftol' : 1e-07}

在 sascha 的代码中。

【讨论】：