Calculating The Inverse of Digamma Function

The digamma function is the logarithmic derivative of the gamma function which is defined for the nonnegative real numbers. When you are working with Beta and Dirichlet distributions, you seen them frequently. Furthermore, if you want to estimate the parameters of a Diricihlet distribution, you need to take the inverse of the digamma function.

I was using the inverse digamma function (invps) of Paul Flacker (code is here). Since my experiments were taking a very long time, I profiled my matlab code, and found out that most of the time the program was taking time in the invps function. I started to search for a better technique and found it on Minka’s technical paper. (This is an excellent report by the way. Thanks.)

Our problem is to find $x$ , given $y$, such that $y=psi(x)$. The solution is
\begin{align}
x^* = \min_x \psi(x)-y
\end{align}
The Newton iteration for this problem is
\begin{align}
x^{(n+1)} = x^{(n)} – \frac{\psi(x^{(n)})-y}{\psi'(x^{(n)})}
\end{align}

In order to minimize the number of Newton iterations, we need to find a good initial value for $x$. The $\psi(x)$ function can be approximated by
\begin{align}
\psi(x) \approx g(x) = \left \{ \begin{matrix} \log(x-1/2) & \text{if } x\geq 0.6 \\ -(1/x+\psi(1)) & \text{if } x<0.6 \end{matrix} \right . \end{align} Therefore we initialize $x$ as \begin{align} x^{(1)} = \left \{ \begin{matrix} \exp(y)+1/2 & \text{if } y\geq -2.22 \\ -1/(y+\psi(1)) & \text{if } y<-2.22 \end{matrix} \right . \end{align}

function Y=invpsi(X)
%INVPSI - Inverse digamma (psi) function.  
%
% The digamma function is the derivative of the log gamma function.  This 
% function calculates the value Y > 0 for a value X such that 
% digamma(Y) = X. The code takes 3 Newton steps after initializing the Y
% according to a good approximation of digamma function.
%
% Source: Thomas Minka, Estimating a Dirichlet distribution, 
%         Tehcnical Report 2012, Appendix C.
%
% Change History :
% Date              Time	   Prog
% 16-Sep-2013      23:40 PM    Baris Kurt	
% Bogazici University, Dept. of Computer Eng. 80815 Bebek Istanbul Turkey
% e-mail : bariskurt@gmail.com    

%initial estimate
M = double(X >= -2.22);
Y = M .*(exp(X) + 0.5) + (1-M) .* -1./(X-psi(1));

% make 3 Newton iterations:
Y = Y - (psi(Y)-X)./psi(1,Y);
Y = Y - (psi(Y)-X)./psi(1,Y);
Y = Y - (psi(Y)-X)./psi(1,Y);

end