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}

### Matlab Code

Here is the matlab implementation. I used only 3 fixed point iterations, which is quite satisfactory.

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

You can also download it as a script invpsi.m