Write Dirichlet Log Likelihood with DCP ruleset

Question

I would like to write the log likelihood of the Dirichlet density as a disciplined convex programming (DCP) optimization problem with respect to the parameters of the Dirichlet distribution alpha. However, the log likelihood

def dirichlet_log_likelihood(p, alpha):
  """Log of Dirichlet density.
  p: Numpy array of shape (K,) that sums to 1.
  alpha: Numpy array of shape (K, ) with positive elements.
  """
  L = np.log(scipy.special.gamma(alpha.sum()))
  L -= np.log(scipy.special.gamma(alpha)).sum()
  L += np.sum((alpha - 1) * np.log(p))
  return L

despite being concave in alpha is not formulated as DCP because it involves the difference of two concave functions np.log(gamma(alpha.sum())) and np.log(gamma(alpha)).sum(). I would like if possible, to formulate this function of alpha so that it follows the DCP ruleset, so that maximum-likelihood estimation of alpha can be performed with cvxpy.

Is this possible, and if so how might I do it?

Answer 1

As you note, np.log(gamma(alpha.sum())) and -np.log(gamma(alpha)).sum() have different curvature, so you need to combine them as

np.log(gamma(alpha.sum()) / gamma(alpha).sum())

to have any chance of modelling them under the DCP ruleset. The combined expression above can be recognized as the logarithm of the multivariate beta function, and since the multivariate beta function can be written as a product of bivariate beta functions (see here), you can expand the log-product to a sum-log expression where each term is of the form

np.log(beta(x,y))

and this is the convex atom you need in your DCP ruleset. What remains of you, to use it in practice, is to feed in an approximation of this atom into cvxpy. The np.log(gamma(x)) approximation here will serve as a good starting point.

Please see math.stackexchange.com for more details.