What is a good optimization algorithm for finding the MAP?

Question

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.

Answer 1

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.

Answer 2

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.