2014-04-23 12 views
7

Tôi đang cố gắng để giảm thiểu các chức năng sau đây với scipy.optimize:giảm thiểu một đa biến, chức năng vi sử dụng scipy.optimize

enter image description here

có độ dốc là thế này:

enter image description here

(cho những người quan tâm, đây là chức năng khả năng của mô hình Bradley-Terry-Luce để so sánh cặp đôi. Rất liên quan chặt chẽ với hồi quy logistic.)

Rõ ràng là việc thêm một hằng số vào tất cả các tham số sẽ không thay đổi giá trị của hàm. Do đó, tôi để cho \ theta_1 = 0. Dưới đây là việc thực hiện các chức năng khách quan và gradient trong python (theta trở thành x đây):

def objective(x): 
    x = np.insert(x, 0, 0.0) 
    tiles = np.tile(x, (len(x), 1)) 
    combs = tiles.T - tiles 
    exps = np.dstack((zeros, combs)) 
    return np.sum(cijs * scipy.misc.logsumexp(exps, axis=2)) 

def gradient(x): 
    zeros = np.zeros(cijs.shape) 
    x = np.insert(x, 0, 0.0) 
    tiles = np.tile(x, (len(x), 1)) 
    combs = tiles - tiles.T 
    one = 1.0/(np.exp(combs) + 1) 
    two = 1.0/(np.exp(combs.T) + 1) 
    mat = (cijs * one) + (cijs.T * two) 
    grad = np.sum(mat, axis=0) 
    return grad[1:] # Don't return the first element 

Dưới đây là một ví dụ về những gì cijs có thể trông giống như:

[[ 0 5 1 4 6] 
[ 4 0 2 2 0] 
[ 6 4 0 9 3] 
[ 6 8 3 0 5] 
[10 7 11 4 0]] 

Đây là mã tôi chạy để thực hiện giảm thiểu:

x0 = numpy.random.random(nb_items - 1) 
# Let's try one algorithm... 
xopt1 = scipy.optimize.fmin_bfgs(objective, x0, fprime=gradient, disp=True) 
# And another one... 
xopt2 = scipy.optimize.fmin_cg(objective, x0, fprime=gradient, disp=True) 

Tuy nhiên, nó luôn luôn thất bại trong phiên đầu tiên:

Warning: Desired error not necessarily achieved due to precision loss. 
     Current function value: 73.290610 
     Iterations: 0 
     Function evaluations: 38 
     Gradient evaluations: 27 

Tôi không thể hiểu tại sao nó không thành công. Lỗi được hiển thị vì dòng này: https://github.com/scipy/scipy/blob/master/scipy/optimize/optimize.py#L853

Vì vậy, "tìm kiếm đường Wolfe" này dường như không thành công, nhưng tôi không biết cách tiến hành từ đây ... Mọi trợ giúp đều được đánh giá cao!

+1

Chức năng gradient của bạn có thể không chính xác. Hãy thử xác minh nó chống lại sự khác biệt hữu hạn (ví dụ: sử dụng [scipy.optimize.check_grad] (http://docs.scipy.org/doc/scipy-dev/reference/generated/scipy.optimize.check_grad.html)) –

+0

@pv. Bạn đặt cược;) Cảm ơn! – lum

Trả lời

2

Như @pv. chỉ ra như một bình luận, tôi đã sai lầm trong việc tính toán gradient. Trước hết, chính xác (toán học) biểu hiện cho gradient của hàm mục tiêu của tôi là:

enter image description here

(chú ý dấu trừ.) Hơn nữa, thực hiện Python của tôi là hoàn toàn sai, vượt quá sai lầm dấu.Dưới đây là dốc cập nhật của tôi:

def gradient(x): 
    nb_comparisons = cijs + cijs.T 
    x = np.insert(x, 0, 0.0) 
    tiles = np.tile(x, (len(x), 1)) 
    combs = tiles - tiles.T 
    probs = 1.0/(np.exp(combs) + 1) 
    mat = (nb_comparisons * probs) - cijs 
    grad = np.sum(mat, axis=1) 
    return grad[1:] # Don't return the first element. 

Để gỡ lỗi nó, tôi đã sử dụng:

  • scipy.optimize.check_grad: cho thấy chức năng dốc của tôi đã được sản xuất kết quả rất xa một xấp xỉ (sai phân hữu hạn) gradient.
  • scipy.optimize.approx_fprime để biết ý tưởng về các giá trị sẽ như thế nào.
  • một vài ví dụ đơn giản được chọn thủ công có thể được phân tích bằng tay nếu cần và một vài truy vấn Wolfram Alpha để kiểm tra sự lành mạnh.
1

Dường như bạn có thể biến đổi nó thành một sự cố ít nhất là vuông. Bằng cách này, bạn phải xác định khoảng thời gian cho mỗi biến số n và số điểm mẫu cho mỗi biến để xây dựng ma trận hệ số.

Trong ví dụ này tôi đang sử dụng cùng số điểm và khoảng thời gian tương tự cho tất cả các biến:

from scipy.optimize import leastsq 
from numpy import exp, linspace, zeros, ones 

n = 4 
npts = 1000 
xs = [linspace(0, 1, npts) for _ in range(n)] 

c = ones(n**2) 

a = zeros((n*npts, n**2)) 
def residual(c): 
    a.fill(0) 
    for i in range(n): 
     for j in range(n): 
      for k in range(npts): 
       a[i+k*n, i*n+j] = 1/(exp(xs[i][k] - xs[j][k]) + 1) 
       a[i+k*n, j*n+i] = 1/(exp(xs[j][k] - xs[i][k]) + 1) 

    return a.dot(c) 

popt, pconv = leastsq(residual, x0=c) 
print(popt.reshape(n, n)) 
#[[ -1.24886411 1.07854552 -2.67212118 1.86334625] 
# [ -7.43330057 2.0935734 37.85989442 1.37005925] 
# [ -3.51761322 -37.49627917 24.90538136 -4.23103535] 
# [ 11.93000731 2.52750715 -14.84822686 1.38834225]] 

EDIT: thêm chi tiết về các hệ số ma trận được xây dựng trên:

enter image description here

+0

Cảm ơn bạn đã cố gắng giúp tôi. Tôi thấy nhiều hơn hoặc ít hơn những gì bạn có nghĩa là, nhưng tôi muốn tránh phù hợp với hình vuông nhỏ nhất. Chức năng mục tiêu của tôi là lồi, vì vậy tôi không thấy lý do gì khiến tôi không thể giảm thiểu nó một cách trực tiếp. – lum

+0

@ lum Tôi thấy điểm của bạn ... dù sao, đây là một giải pháp rất mạnh mẽ trong trường hợp bạn cần nó .. –

Các vấn đề liên quan