Método para encontrar una solución estrictamente mayor que cero para la Progtwigción lineal Scipy de Python

Scipy NNLS realiza esto:

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

¿Cuál es la forma alternativa de hacerlo si busco una solución estrictamente distinta de cero (es decir, x > 0 )?

Aquí está mi código LP usando el NNLS de Scipy:

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

El uso:

 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 la solución cero arriba.

Debería preguntar si realmente quiere x > 0 lugar de x >= 0 . Por lo general, la última restricción sirve para dispersar el resultado, y los ceros en x son deseables. Aparte de eso, las restricciones son prácticamente equivalentes.

Si restringe x para que sea estrictamente mayor que cero, los 0 simplemente se convertirán en números positivos muy pequeños. Si la alma se pudiera mejorar con valores más grandes, también obtendría estos valores con la restricción original.

Demuestre esto definiendo la siguiente optimización: Resuelva argmin_x || Ax - b ||_2 for x>=eps argmin_x || Ax - b ||_2 for x>=eps . Mientras que eps > 0 esto también satisface x > 0 . Mirando la x resultante para diferentes eps , obtenemos:

resultado

Lo que ves es que para mall eps casi no hay diferencia en la función objective y x[1] (uno de los 0 en la solución original) se acerca cada vez más a 0. Así, el paso infinitesimal de x>0 a x>=0 apenas cambia nada en la solución. A efectos prácticos son totalmente similares. Sin embargo, x>=0 tiene la ventaja de que obtiene 0s reales en lugar de 1.234e-20, lo que ayuda a dispersar la solución.

Aquí está el código para la ttwig anterior:

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

PS He intentado implementar la restricción x > 0 en mi código con lambda x: [1 if i>0 else -1 for i in x] , pero no puede converger.

Si está realmente seguro de que desea soluciones estrictamente positivas, puede usar lsq_linear, que está disponible en la versión scipy más reciente. Permite un poco más de control sobre los límites que 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])