为Python的Scipy线性编程找到严格大于零的解决方案的方法 [英] Method for finding strictly greater-than-zero solution for Python's Scipy Linear Programing
问题描述
Scipy NNLS 执行以下操作:
Solve argmin_x || Ax - b ||_2 for x>=0.
如果我寻求帮助，还有什么替代方法?
严格非零的解决方案(即x > 0)?
What's the alternative way to do it if I seek
strictly non-zero solution (i.e. x > 0) ?
这是我使用Scipy的NNLS的LP代码:
Here is my LP code using Scipy's NNLS:
import numpy as np
from numpy import array
from scipy.optimize import nnls
def by_nnls(A=None, B=None):
""" Linear programming by NNLS """
#print "NOF row = ", A.shape[0]
A = np.nan_to_num(A)
B = np.nan_to_num(B)
x, rnorm = nnls(A,B)
x = x / x.sum()
# print repr(x)
return x
B1 = array([ 22.133, 197.087, 84.344, 1.466, 3.974, 0.435,
8.291, 45.059, 5.755, 0.519, 0. , 30.272,
24.92 , 10.095])
A1 = array([[ 46.35, 80.58, 48.8 , 80.31, 489.01, 40.98,
29.98, 44.3 , 5882.96],
[ 2540.73, 49.53, 26.78, 30.49, 48.51, 20.88,
19.92, 21.05, 19.39],
[ 2540.73, 49.53, 26.78, 30.49, 48.51, 20.88,
19.92, 21.05, 19.39],
[ 30.95, 1482.24, 100.48, 35.98, 35.1 , 38.65,
31.57, 87.38, 33.39],
[ 30.95, 1482.24, 100.48, 35.98, 35.1 , 38.65,
31.57, 87.38, 33.39],
[ 30.95, 1482.24, 100.48, 35.98, 35.1 , 38.65,
31.57, 87.38, 33.39],
[ 15.99, 223.27, 655.79, 1978.2 , 18.21, 20.51,
19. , 16.19, 15.91],
[ 15.99, 223.27, 655.79, 1978.2 , 18.21, 20.51,
19. , 16.19, 15.91],
[ 16.49, 20.56, 19.08, 18.65, 4568.97, 20.7 ,
17.4 , 17.62, 25.51],
[ 33.84, 26.58, 18.69, 40.88, 19.17, 5247.84,
29.39, 25.55, 18.9 ],
[ 42.66, 83.59, 99.58, 52.11, 46.84, 64.93,
43.8 , 7610.12, 47.13],
[ 42.66, 83.59, 99.58, 52.11, 46.84, 64.93,
43.8 , 7610.12, 47.13],
[ 41.63, 204.32, 4170.37, 86.95, 49.92, 87.15,
51.88, 45.38, 42.89],
[ 81.34, 60.16, 357.92, 43.48, 36.92, 39.13,
1772.07, 68.43, 38.07]])
用法:
In [9]: by_nnls(A=A1,B=B1)
Out[9]:
array([ 0.70089761, 0. , 0.06481495, 0.14325696, 0.01218972,
0. , 0.02125942, 0.01906576, 0.03851557]
请注意上面的零解.
推荐答案
您应该质疑您是否真的想要x > 0而不是x >= 0.通常，后一个约束用于稀疏结果，并且x中的零是可取的.除此之外，约束条件实际上是等效的.
You should question if you really want x > 0 instead of x >= 0. Usually the latter constraint serves to sparsify the result, and zeros in x are desirable. Apart from that, the constraints are practically equivalent.
如果将x严格限制为大于零，则0将仅变为非常小的正数.如果可以通过更大的值来改善解决方案，那么您也将获得具有原始约束的这些值.
If you constrain x to be strictly greater than zero, the 0s will simply become very very small positive numbers. If the soulution could be improved by larger values, you would get these values with the original constraint too.
让我们通过定义以下优化对此进行演示:解决argmin_x || Ax - b ||_2 for x>=eps.当eps > 0时，它也满足x > 0.查看不同的eps的结果x，我们得到:
Let's demonstrate this by defining the following optimization: Solve argmin_x || Ax - b ||_2 for x>=eps. While eps > 0 this also satisfies x > 0. Looking at the resulting x for different eps, we get:
您看到的是对于购物中心eps，目标函数几乎没有差异，x[1](原始解决方案中的0之一)越来越接近0.
因此，从x>0到x>=0的无穷小步几乎不会改变解中的任何内容.出于实际目的，它们是完全相似的.但是，x>=0的优点是您得到的是实际0，而不是1.234e-20，这有助于简化解决方案.
What you see is that for mall eps there is hardly any difference in the objective function and x[1] (one of the 0s in the original solution) gets closer and closer to 0.
Thus, the infinitesimal step from x>0 to x>=0 hardly changes anything in the solution. For practical purposes they are totally similar. However, x>=0 has the advantage that you get actual 0s rather than 1.234e-20, which helps in sparsifying the solution.
这是上面剧情的代码:
from scipy.optimize import fmin_cobyla
import matplotlib.pyplot as plt
def by_minimize(A, B, eps=1e-6):
A = np.nan_to_num(A)
B = np.nan_to_num(B)
def objective(x, A=A, B=B):
return np.sum((np.dot(A, x) - B)**2)
x0 = np.zeros(A.shape[1])
x = fmin_cobyla(objective, x0, lambda x: x-eps)
return x / np.sum(x), objective(x)
results = []
obj = []
e = []
for eps in np.logspace(-1, -6, 100):
x, o = by_minimize(A=A1, B=B1, eps=eps)
e.append(eps)
results.append(x[1])
obj.append(o)
h1 = plt.semilogx(e, results, 'b')
plt.ylabel('x[1]', color='b')
plt.xlabel('eps')
plt.twinx()
h2 = plt.semilogx(e, obj, 'r')
plt.ylabel('objective', color='r')
plt.yticks([])
P.S.我试图用lambda x: [1 if i>0 else -1 for i in x]在我的代码中实现x > 0约束，但是无法收敛.
P.S. I have tried to implement the x > 0 constraint in my code with lambda x: [1 if i>0 else -1 for i in x], but it fails to converge.
这篇关于为Python的Scipy线性编程找到严格大于零的解决方案的方法的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
