Peace be upon you,
At the end of this question, I have shown that how computing MLE on an i.i.d Beta distributed data, results in the following system \begin{align*} &\begin{cases} \psi(\alpha)-\psi(\alpha+\beta) = c_1\\ \psi(\beta)-\psi(\alpha+\beta) = c_2 \end{cases}\\ &c_1 = \frac{1}{n} \left[ \sum_{i=1}^n {log(x_i)}\right] &c_2 = \frac{1}{n} \left[ \sum_{i=1}^n {log(1-x_i)} \right] \end{align*} where $\psi(.)$ is the digamma function and $x_i\in(0,1]$.
Since I have about 19k of the above system to be solved, then using the iterative methods like Newton-Raphson would take a long time. So, I would like to find a fast-computable late-diverging approximating answer for it. I followed the idea of estimating the $\psi(.)$ function by a, likely, invertible function and substituting it and then solving a new system which results in a closed-form approximating solution.
The simplest candidate for this idea is the Taylor expansion; but it would have a very bad accuracy by getting far from the origin or in other words it diverges quickly.
I suppose that for some functions the divergence would occur very late. For example, $\ln(x)$ is very similar to $ψ(x)$ and it may keep the precision in a large domain. but for such suggestion there is another problem: how to solve the system after substitution.
Maybe, someone will suggest another way for solving the above system. But if some user wants to solve it by the aforementioned technique, he/she can find a much similar basis (than $\{1,x,x^2,...\}$) for a $ψ(x)$ estimation which makes the above system has a closed-form answer.
Note that I need a "late-diverging" fast-computable answer, either a closed-form approximating answer or any other fast-computable one.
But, how did I reach this system?
We know that the probability density function of the Beta distribution is \begin{align*} f(x|\alpha, \beta)=\frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha,\beta)} \end{align*} So, the likelihood is \begin{align*} \mathcal{L}(\alpha,\beta)=f(x_1,x_2,...,x_n|\alpha,\beta)=\prod_{i=1}^n{f(x_i|\alpha,\beta)} \end{align*} For pre-simplifying the computations, we proceed to find the maximum of the log-likelihood \begin{align*} log(\mathcal{L}(\alpha,\beta))=\sum_{i=1}^n{log\left(f(x_i|\alpha,\beta)\right)} \end{align*} For maximizing the above log-likelihood, we should solve the following system \begin{align*} \begin{cases} \frac{\partial}{\partial \alpha}\left[ \sum_{i=1}^n {(\alpha-1)log(x_i)+(\beta-1)log(1-x_i)-log(B(\alpha,\beta))} \right]=0\\ \frac{\partial}{\partial \alpha}\left[ \sum_{i=1}^n {(\alpha-1)log(x_i)+(\beta-1)log(1-x_i)-log(B(\alpha,\beta))} \right]=0 \end{cases} \end{align*} So, we can write \begin{align*} \begin{cases} \left[ \sum_{i=1}^n {log(x_i)}\right] =n\frac{\partial}{\partial\alpha}\left[log(B(\alpha,\beta)\right]\\ \left[ \sum_{i=1}^n {log(1-x_i)} \right]=n\frac{\partial}{\partial\beta}\left[log(B(\alpha,\beta)\right] \end{cases} \end{align*} Using the formula for partial derivation of Beta function we can write \begin{align*} \begin{cases} \left[ \sum_{i=1}^n {log(x_i)}\right] =n\left[\psi(\alpha)-\psi(\alpha+\beta)\right]\\ \left[ \sum_{i=1}^n {log(1-x_i)} \right]=n\left[\psi(\beta)-\psi(\alpha+\beta)\right] \end{cases} \end{align*} where $\psi(.)$ is the digamma function. we can rewrite the above system as \begin{align*} &\begin{cases} \psi(\alpha)-\psi(\alpha+\beta) = c_1\\ \psi(\beta)-\psi(\alpha+\beta) = c_2 \end{cases}\\ &c_1 = \frac{1}{n} \left[ \sum_{i=1}^n {log(x_i)}\right] &c_2 = \frac{1}{n} \left[ \sum_{i=1}^n {log(1-x_i)} \right] \end{align*}