The Davidon-Fletcher-Powell method has worked well for me in problems like this.
Below is some super hacky code that I wrote for using this with a line-search method for an Econometrics class. You'll have to implement your own derivatives and log-likelihood wrapper functions to get it to work, but maybe it will be helpful. It should only depend on NumPy and SciPy.
def line_search_iteration(data,gamma,pi,v_vector,c,A):
# Calculate the derivative at the current value of beta
diff = -1.0*compute_numerical_derivative(data,gamma,pi,v_vector,c)
# Update normalized direction
d = -A.dot(diff)
d_n = 1.0/(1+np.linalg.norm(d))*d
# Compute a bracket of the root.
# Compute function values to determine a root bracket.
lam_l = 0.0
lam_h = 1.0
Qh = []
tmp0 = np.asarray([[gamma[0,0]],[gamma[1,0]],[gamma[2,0]],[gamma[3,0]],[pi]])
tmpL = tmp0 + lam_l*d_n
tmpH = tmp0 + lam_h*d_n
Qh.append(-compute_log_likelihood(data,tmpL[0:-1,:],tmpL[-1,0],v_vector) )
Qh.append(-compute_log_likelihood(data,tmpH[0:-1,:],tmpH[-1,0],v_vector) )
# Search for a root bracket.
while Qh[-1] >= Qh[-2]:
lam_l = lam_l + (lam_h-lam_l)
lam_h = 2.0*lam_h
tmpH = tmp0 + lam_h*d_n
Qh.append(-compute_log_likelihood(data,tmpH[0:-1,:],tmpH[-1,0],v_vector) )
# Use golden search to compute optimal length
# Do golden section to obtain the optimal value of lam
# 0th iteration
lam_m1 = lam_l + 0.382*(lam_h - lam_l)
lam_m2 = lam_l + 0.618*(lam_h - lam_l)
tmp_m1 = tmp0 + lam_m1*d_n
tmp_m2 = tmp0 + lam_m2*d_n
Q_m1 = -compute_log_likelihood(data,tmp_m1[0:-1,:],tmp_m1[-1,0],v_vector)
Q_m2 = -compute_log_likelihood(data,tmp_m2[0:-1,:],tmp_m1[-1,0],v_vector)
while (lam_h-lam_l) >= 0.0001:
if Q_m2 > Q_m1:
# lam_l doesn't change
lam_h = np.copy(lam_m2)
lam_m2 = np.copy(lam_m1)
lam_m1 = lam_l + 0.382*(lam_h-lam_l)
tmp_m1 = tmp0 + lam_m1*d_n
tmp_m2 = tmp0 + lam_m2*d_n
Q_m1 = -compute_log_likelihood(data,tmp_m1[0:-1,:],tmp_m1[-1,0],v_vector)
Q_m2 = -compute_log_likelihood(data,tmp_m2[0:-1,:],tmp_m1[-1,0],v_vector)
lam_avg = 0.5*(lam_l + lam_h)
else:
# lam_h doesn't change
lam_l = np.copy(lam_m1)
lam_m1 = np.copy(lam_m2)
lam_m2 = lam_l + 0.618*(lam_h-lam_l)
tmp_m1 = tmp0 + lam_m1*d_n
tmp_m2 = tmp0 + lam_m2*d_n
Q_m1 = -compute_log_likelihood(data,tmp_m1[0:-1,:],tmp_m1[-1,0],v_vector)
Q_m2 = -compute_log_likelihood(data,tmp_m2[0:-1,:],tmp_m1[-1,0],v_vector)
lam_avg = 0.5*(lam_l + lam_h)
# Now, return the necessary components for making the update.
return diff, d_n, lam_avg
####################
""" Function for carrying out a single DFP iteration.
"""
def dfp_iteration(data,gamma,pi,v_vector,c,A):
# Get the necessary step direction and optimal step length for the derivative
diff, d_n, lam = line_search_iteration(data,gamma,pi,v_vector,c,A)
# Update the estimate of beta and the scaling matrix A.
p = lam*d_n
gamma_aug = np.asarray([[gamma[0,0]],[gamma[1,0]],[gamma[2,0]],[gamma[3,0]],[pi]])
gamma_aug = gamma_aug + p
err_prime = np.linalg.norm(p)
gamma_new = gamma_aug[0:-1,:]
pi_new = gamma_aug[-1,0]
q = -1*compute_numerical_derivative(data,gamma_new,pi_new,v_vector,c) - diff
num1 = p.dot(np.transpose(p))
den1 = np.transpose(p).dot(q)
num2 = A.dot(q.dot(np.transpose(q).dot(A)))
den2 = np.transpose(q).dot(A.dot(q))
A = A + num1/den1 - num2/den2
err = np.transpose(diff).dot(diff)
return gamma_new, pi_new, A, err, err_prime
There is also CVXOpt for Python, but I've never used it. Not sure how adaptable it is for your problem.