Method for finding strictly greater-than-zero solu

问题:

Scipy NNLS perform this:

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

What's the alternative way to do it if I seek strictly non-zero solution (i.e. x > 0) ?

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

The usage:

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]

Note the zero solution above.

回答1:

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.

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.

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:

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.

Here is the code for above plot:

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. 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.

回答2:

If you are really sure you want strictly positive solutions, you can use lsq_linear which is available in the most recent scipy version. It allows a bit more control over the bounds than 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])