Currently, I encountered a problem of approximating the following
series:


$$
(I-X)^{-\frac{1}{2}}=I+\frac{1}{2}X+\frac{1\cdot3}{2\cdot4}X^{2}+\frac{1\cdot3\cdot5}{2\cdot4\cdot6}X^{3}+\ldots
$$


where $X$ is a diagonalizable matrix and the largest (in absolute
value) eigenvalue is less than $1-\frac{1}{\kappa},\kappa>1$. What
I am looking for is something like the following


$$
(1-X)^{-1}=I+X+X^{2}+X^{3}+\dots\approx\prod\_{k=0}^{d}(1+X^{2^{k}})
$$


With $O(\log(\kappa\log(1/\epsilon)))$ multiplications and summations,
we can achieve accuracy of $O(\epsilon)$.


We can assume the following representations are cheap to get (we assume
that their computational cost is $1$):


1. $1+aX^{n},a\in(-1,1),n\in[-1,0,1,\dots)$
2. Multiplication between any two cheap representations.
3. Summation between any two cheap representations.


So how can we approximate $(I-X)^{-\frac{1}{2}}$ efficiently?


My initial idea is 


$$
(I-X)^{-\frac{1}{2}}=I+\frac{1}{2}X(1+\frac{3}{4}X(1+\dots(1+\frac{2n-1}{2n}X(1-X)^{-1}))\dots)
$$


But it is not efficient enough.



Approximate the square root of (1-X) efficiently through (nested) products
Not exactly the same thing, but you can use Newton's method easily in your framework. As you can see for instance [here](https://en.wikipedia.org/wiki/Fast_inverse_square_root#Newton.27s_method), Newton's method for the inverse square root can be computed with **additions and multiplications only**, and the same algebraic expression that works for scalars can be used for diagonalizable matrices without modifications (proof: in the eigenvector basis, doing the iteration on the matrix level is equivalent to running $n$ scalar iterations inside a diagonal matrix).


It converges quadratically, so it should have the convergence rate that you need.


By the way, your initial idea is known in numerical analysis as [Horner's scheme](https://en.wikipedia.org/wiki/Horner_method) for polynomial evaluation.


