Python的Scipy线性规划寻找严格大于零解的方法答案

【问题标题】：Method for finding strictly greater-than-zero solution for Python's Scipy Linear ProgramingPython的Scipy线性规划寻找严格大于零解的方法
【发布时间】：2025-11-23 11:20:04
【问题描述】：

Scipy NNLS 执行此操作：

Solve argmin_x || Ax - b ||_2 for x>=0.

如果我寻求另一种方法来做到这一点严格非零解（即x > 0）？

这是我使用 Scipy 的 NNLS 的 LP 代码：

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]

注意上面的零解。

【问题讨论】：

为什么要严格大于零的系数？如果您试图施加这样的约束，那么它们可能会被推向任意小的正数。值 0 和 2⁻¹⁰²²（这是 64 位浮点数的最小可用非零值）之间是否有任何有意义的差异？

标签： python numpy scipy linear-algebra mathematical-optimization

【解决方案1】：

如果你真的确定你想要严格积极的解决方案，你可以使用最新的 scipy 版本中提供的lsq_linear。与nnls 相比，它允许对边界进行更多控制。

In [37]: from scipy.optimize import lsq_linear

In [38]: lsq_linear(A1, B1, bounds=(0.001, np.inf))
Out[38]: 
 active_mask: array([ 0, -1,  0,  0, -1, -1,  0,  0,  0])
        cost: 3784.3150152135881
         fun: array([ -0.06189388, -56.45892624,  56.28407376,   2.97647016,
         0.46847016,   4.00747016,  18.24947887, -18.51852113,
         0.19599207,   7.32663679,  15.0829264 , -15.1890736 ,
        -0.14570891,  -0.24341795])
     message: 'The first-order optimality measure is less than `tol`.'
         nit: 17
  optimality: 5.4491449547056092e-11
      status: 1
     success: True
           x: array([ 0.05506904,  0.001     ,  0.00501077,  0.01112669,  0.001     ,
        0.001     ,  0.00154812,  0.00147833,  0.00300156])

【讨论】：

【解决方案2】：

您应该质疑您是否真的想要x > 0 而不是x >= 0。通常后一个约束用于稀疏化结果，并且 x 中的零是可取的。除此之外，约束实际上是等效的。

如果您将 x 限制为严格大于零，则 0 将简单地变成非常非常小的正数。如果可以通过更大的值来改进解决方案，那么您也可以在原始约束条件下获得这些值。

让我们通过定义以下优化来证明这一点：求解argmin_x || Ax - b ||_2 for x>=eps。而eps > 0 这也满足x > 0。查看不同eps 的结果x，我们得到：

您看到的是，对于 mall eps，目标函数几乎没有任何区别，而 x[1]（原始解决方案中的 0 之一）越来越接近 0。因此，从x>0 到x>=0 的微小步骤几乎不会改变解决方案中的任何内容。出于实际目的，它们完全相似。但是，x>=0 的优势在于您得到的是实际的 0 而不是 1.234e-20，这有助于解决方案的稀疏化。

这是上图的代码：

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([])

附：我尝试使用lambda x: [1 if i>0 else -1 for i in x] 在我的代码中实现x > 0 约束，但它未能收敛。

【讨论】：