# What is a good optimization algorithm for finding the MAP?

I'm working on a Bayesian inference package and I have a function that finds the maximum a posteriori point. Right now I'm using scipy's implementation of BFGS to find that minimum of the log posterior, but I've had some trouble with handling large negatives and -infinities. I'd like to see if there's a better algorithm, particularly something that works well on problems that

1. Will be stiff in some dimensions (like scale parameters)
2. Have boundaries (like uniform distributions)
3. Have a largish number of dimensions (say at least 100)

I'm looking for suggestions on what optimization algorithms to look at.

I'm okay with implementing something myself.

Three suggestions:

1. Would you consider looking CERES solver by Google?

2. Being somewhat naive on this, would you consider using a Stochastic Gradient descent, possibly with some form of regularization (like the one done here)?

3. Check NLopt and in particular the implementations of COBYLA and BOBYQA (and possibly NEWUOA). BOBYQA practically forms quadratic models by interpolation and allows for box constraints (bounds) on the parameters. (COBYLA uses Linear Approximations)

Also, just because you mentioned -Inf(s); I suspect you get those out of the "determinant" term. How exactly are you computing the determinants? For example, just using log(det(A)) in Matlab gave me -Inf(s) in some cases that 2*sum(log(diag(chol(A)))) actually gave me a "number". Maybe you would like to check the numerical stability of the matrix decomposition routines you employ.

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.