What is the most efficient way to solve system of equations containing the digamma function?

Question

What is the most efficient way to solve system of equations involving the digamma function?

I have a vector v and I want to solve for a vector w such that for all i:

digamma(sum(w)) - digamma(w_i) = v_i

and

w_i > 0

I found the gsl function gsl_sf_psi, which is the digamma function (calculated using some kind of series.) Is there an identity I can use to reduce the equations? Is my best bet to use a solver? I am using C++0x; which solver is easiest to use and fast?

---

From my preliminary research, digamma is not easily invertible (searching for inverse digamma gives algorithms that work by binary search), so it makes sense that there would be no simplification to the whole system.

So, using a solver now leaves two problems: dealing with the fact that digamma is very slow to calculate, and dealing with the restriction that w_i > 0, or else digamma(w_i) will crash for w_i = 0.

For the first problem, I thought maybe I should implement a cache for recently calculated values of digamma -- I thought that would be a good idea, but don't know much about how root-finding algorithms work.

My idea was to solve the second problem was to find w'_i = log(w_i). That way, w'_i are on the whole line. I wonder if this is a good idea. There's probably no function to find digamma(exp(w')) directly? Also, the algorithm might take steps in w' space and not improve things because the mapping from w'->w loses some precision and so two elements of w' might map to the same w.

There's still the question of finding a good, fast rootfinding algorithm. I think I may ask that in a separate question.

Thanks...

Answer 1

I would suggest that using a solver would be the best idea, mainly because considering the various stability and convergence regions for various equations might be tricky and it is no use reinventing the wheel. While I have never really solved a system like the one you mentioned, I think one of the following libraries are highly likely to have a solution that you want:

• http://www.gnu.org/software/gsl/
• http://www.scimath.com/
• http://www.moshier.net/#Cephes
• There are more that you can google for certainly :) [These are what I a familiar to]

Also, if neither of these has exactly what you want, you could see how GNU Octave or something alike solve the system and then read their documentation about the algorithm they use to implement the functions needed to solve it. From there it is more about figuring out how the algorithm is used, how to implement it, and in which cases is it worthwhile (documentation of Octave, Matlab, Mathematica is very comprehensive and lists the publications that have the algorithm defined in most cases), there are also Scilab and SageMath if you are looking for opensource/free alternatives, and there are ways to use routines from these from within C++ (but I'm not sure how easy or difficult that would be)

Hope that helps.