Cvxpy: Missing atom: gammaln

[Bonnevie]

I heard Boyd joking at a recent talk that you never get requests for new atoms, so here is one:
I am a bit hamstringed by the lack of gammaln (logarithmic gamma function) in your list of atoms. It's definitely convex (e.g. see https://proofwiki.org/wiki/Log_of_Gamma_Function_is_Convex_on_Positive_Reals or just plot the second derivative), so I think it would be a fruitful addition to cvxpy. It occurs with some frequency when working with probabilistic models as many exponential family distributions have a gamma function embedded in their partition function.
By extension, the multivariate gamma (multigammaln in scipy) would also be useful.

[SteveDiamond]

I'd suggest approximating the log_gamma function as a piecewise linear function. If you find a really good approximation we can add it to cvxpy. This is how log_norm_cdf works in CVX.

[Bonnevie]

I am told the Lanczos approximation is the way to go with gammaln. Do you need an implementation in C? It's available in the cmath package as "lgamma".

[SteveDiamond]

The point of approximating the function isn't to evaluate it. It's to create a representation that CVXPY can work with and incorporate into the final form of the problem. Specifically, lgamma needs to be represented as partial optimization over a cone program, e.g.,

lgamma(x) = minimize_{over t} c^Tt
subject to dx + At + b in K

where t in R^n is a new variable, d in R^m, A in R^{m x n} and b in R^m are known constants, and K is a convex cone of a supported type (nonnegative and equal to 0 are the easiest to think about).

[SteveDiamond]

I'm closing this. Feel free to reopen it if there's more interest/progress that you want to contribute.

[cdipaolo]

Hi @Bonnevie @SteveDiamond, I also ran into this issue (in the same problem space regarding exponential family ML estimates) and wanted to share the approximations I came up with to get around it. I'm sure there's something much more sophisticated that would work better, but this appears to work well enough and might help someone who runs into this issue later. This doesn't feel robust enough to add directly in to the library, but perhaps a future iteration of it would work for that purpose.

The short of it is that you can use the following expression as a drop-in log-gamma replacement in CVXPY:

x = cp.Variable()
approx_log_gamma = lambda x: cp.maximum(
    2.18382 - 3.62887*x,
    1.79241 - 2.4902*x,
    1.21628 - 1.37035*x,
    0.261474 - 0.28904*x,
    0.577216 - 0.577216*x,
    -0.175517 + 0.03649*x,
    -1.27572 + 0.621514*x,
    -0.845568 + 0.422784*x,
    -0.577216*x - cp.log(x),
    0.918939 - x - cp.entr(x) - 0.5*cp.log(x),
)
approx_log_gamma(x) # Expression(CONVEX, UNKNOWN, ())

This function was constructed by taking maximum of the asymptotic expansion at zero and the asymptotic expansion at infinity, as well as the tangent planes at the two zeros (1 and 2) in order to uniformly bound the relative and absolute error in the approximation. At this point, I greedily added in tangent planes at places where the absolute/relative error was maximized until the absolute error was under 0.05 uniformly and the relative error felt small enough. As a result of this, the approximation lower bounds the true log-gamma function uniformly, and has zero error as you approach 0 or infinity.

The maximum relative error in approximation of the above function is 0.266, while the maximum absolute error is 0.036. See attached plot and notebook for more details there.

To test how this works in practice for this problem domain, I ran a simple experiment where I computed the maximum likelihood estimate of the shape parameter of a gamma distribution (which involves a log-gamma term in the objective). For each of 40 shape and 40 scale parameters non-uniformly taken from the range [0.1,100], I sampled 100 values from the resulting gamma distribution, and computed the maximum likelihood estimate with scale fixed using CVXPY, as well as Newton iteration with the initial guess for the Newton iteration given as in the wiki page [1].

|Approximation Used|75-th Perc. Relative Error| Max Relative Error| Mean Relative Error| Median Relative Error|
|----|------|-----|-----|-----|
|Just Asymptotics | 0.011320 | 0.139755 | 0.014509 | 0.00007|
|Full Approximation|0.006321|0.190309|0.024013|0.00007|

The "Just Asymptotics" approximation drops the tangent plane terms in the approximation above and just uses the maximum of the nonlinear expansions at 0 and infinity, while "Full Approximation" is exactly the one presented above. For these situations where statistical estimation error likely well-outweighs the approximation error above, this seems to be a good enough, and depending on the need you might want to use the simpler approximation or the more complex one. Of course, in high-precision areas both of these might not suffice.

Notebook w/ Code to Reproduce Results: Approximating the Log-Gamma Function.ipynb.zip

[SteveDiamond]

This is fantastic! Thanks so much for doing this, and sharing it with us! I'll let you know if we add your approximation. I'm definitely open to it.