2015-10-28 16 views
5

非負最小二乗を解くためのコードは次のとおりです。 scipy.nnls.Scipy NNLS関数解に制約を含める方法1

import numpy as np 
from scipy.optimize import nnls 

A = np.array([[60, 90, 120], 
       [30, 120, 90]]) 

b = np.array([67.5, 60]) 

x, rnorm = nnls(A,b) 

print x 
#[ 0.   0.17857143 0.42857143] 
# Now need to have this array sum to 1. 

が1にそれの和が、私はそれをどのように行うことができるようにxソリューションに制約を適用するために私は何をしたいです使用 ?

答えて

6

nnlsを直接使用してFortran codeという名前で呼び出すことはできません。余分な制約はありません。しかし、いずれかの式の和が第三の方程式として導入することができるという制約が、これは今一次方程式の集合であるようにあなたの例のシステムは、フォームの

60 x1 + 90 x2 + 120 x3 = 67.5 
30 x1 + 120 x2 + 90 x3 = 60 
    x1 +  x2 +  x3 = 1 

あり、厳密解を求めることができますx=np.dot(np.linalg.inv(A),b)からx=[0.6875, 0.3750, -0.0625]となります。これには、が負である必要があります。したがって、xがこの問題に肯定的である場合、正確な解は存在しません。 xが陽性であることが拘束され、これは使用して得られる近似解について

import numpy as np 
from scipy.optimize import nnls 

#Define minimisation function 
def fn(x, A, b): 
    return np.sum(A*x,1) - b 

#Define problem 
A = np.array([[60., 90., 120.], 
       [30., 120., 90.], 
       [1., 1., 1. ]]) 

b = np.array([67.5, 60., 1.]) 

x, rnorm = nnls(A,b) 

print(x,x.sum(),fn(x,A,b)) 

x.sum()=0.95x=[0.60003332, 0.34998889, 0.]を与えます。

私はあなたが和制約を含む、より一般的な解決策を望んでいた場合、あなたは以下の形式で明示的な制約/限界と最小限使用する必要があると思う、x=[0.674999366, 0.325000634, 0.]x.sum()=1与え

import numpy as np 
from scipy.optimize import minimize 
from scipy.optimize import nnls 

#Define problem 
A = np.array([[60, 90, 120], 
       [30, 120, 90]]) 

b = np.array([67.5, 60]) 

#Use nnls to get initial guess 
x0, rnorm = nnls(A,b) 

#Define minimisation function 
def fn(x, A, b): 
    return np.linalg.norm(A.dot(x) - b) 

#Define constraints and bounds 
cons = {'type': 'eq', 'fun': lambda x: np.sum(x)-1} 
bounds = [[0., None],[0., None],[0., None]] 

#Call minimisation subject to these values 
minout = minimize(fn, x0, args=(A, b), method='SLSQP',bounds=bounds,constraints=cons) 
x = minout.x 

print(x,x.sum(),fn(x,A,b)) 

。最小から、合計は正しいですが、xの値はnp.dot(A,x)=[ 69.75001902, 59.25005706]では正しくありません。

+0

あなたの最初のコードブロック 'print(x、x.sum()、fn(x、A、b))' 'fn'はどこから来ますか? – neversaint

+1

申し訳ありませんが、2番目の例(Ax-b)のfnである必要があります。私は訂正しました。 –

+0

ありがとうございます。私の他の関連する[質問]を見てケア(http://stackoverflow.com/questions/33404908/best-way-to-scale-the-matrix-variables-in-scipy-linear-programming-scheme)? – neversaint

関連する問題