Form Dirichlet Log Density with DCP ruleset.

Question

I would like to formulate the log likelihood of the Dirichlet density as a disciplined convex programming (DCP) optimization problem with respect to $\alpha \in \mathbb{R}^K_{\succ 0}$. The log likelihood

$$\mathcal{L}(\textbf{p}, \alpha) = \log \frac{\Gamma \left(\sum_k \alpha_k \right)}{\prod_k \Gamma (\alpha_k)} \prod_k p_k^{\alpha_k - 1}$$ $$= \log \Gamma \left(\sum_k \alpha_k \right) - \sum_k \log \Gamma(\alpha_k) + \sum_k (\alpha_k - 1) \log p_k $$

despite being concave in $\alpha$ is not formulated as DCP because is involves the difference of $\log \Gamma \left(\sum_k \alpha_k \right)$ and $\sum_k \log \Gamma(\alpha_k)$ which each have positive curvature. Note that despite this, the function is concave because $\log \Gamma (\alpha)$ is convex on $\alpha \succ 0$ and thus

$$\log \Gamma \left(\sum_k \alpha_k \right) \leq \sum_k \log \Gamma(\alpha_k)$$

holds from Jensen's inequality, which is sufficient to show that their difference is concave. However, in order to use certain convex optimization solvers such as cvxpy, I need to formulate this function using the DCP ruleset, which does not allow the difference of two functions with the same curvature in general. Is it possible to do so, and what approaches might I take to achieve this?

Answer 1

Let $$ \beta(x_1,x_2) = \frac{\Gamma \left( x_1 \right) \Gamma \left( x_2 \right)}{\Gamma \left(x_1 + x_2 \right)} $$

be the beta function, and

$$ B(x) = \frac{\prod_{j=1}^n \Gamma \left( x_j \right)}{\Gamma \left(\sum_{j=1}^n x_j \right)}$$

be the multivariate beta function, and note their beautiful relation (just expand and see how everything cancels out to give you exactly what you expect):

$$ B(x) = \prod_{k=2}^{n} \beta \left(\sum_{j=1}^{k-1} x_j,x_k \right) $$

What you are asking for is a way to represent the concave expression $$\log \left( B(x)^{-1} \prod_{j=1}^n p_j^{(x_j - 1)} \right)$$

which, since $\log(x^{-1}) = -\log(x)$, is equivalent to representing the convex expression: $$\log \left( B(x) \prod_{j=1}^n p_j^{(1 - x_j)} \right)\\ = \sum_{k=2}^{n} \log \left( \beta \left(\sum_{j=1}^{k-1} x_j,x_k \right) \right) + \sum_{j=1}^n \log \left(p_j^{(1 - x_j)} \right)\\ = \sum_{k=2}^{n} \log \left( \beta \left(\sum_{j=1}^{k-1} x_j,x_k \right) \right) + \sum_{j=1}^n \log (p_j) (1 - x_j)\\ = \sum_{k=2}^{n} y_k + \sum_{j=1}^n \log (p_j) (1 - x_j) $$

where $$ y_k \geq \log \left( \beta \left(\sum_{j=1}^{k-1} x_j,x_k \right) \right) \quad \mbox{ for } k=2,\ldots,n. $$

In summary, the DCP representation costs one linear expression plus $n-1$ convex inequalities representable by this atom: $$\log \left( \beta \left( x_1,x_2 \right) \right)$$

This atom is convex by logarithmic convexity of the bivariate beta function as shown, e.g., in this paper:

• D. S. Dragomir, R.P. Agarwal, and N. S. Barnett, Inequalities for Beta and Gamma functions via some classical and new integral inequalities, Journal of Inequalities and Applications,5(2000), no. 2, 103–165.

Unfortunately, this atom is neither available in cvxpy nor treated in the MOSEK cookbook, so what remains of you is to device an approximation of it. You can find inspiration, e.g., in the approximation for $\log \left( \Gamma (x) \right)$ found in the cvxpy issue tracker.