The matrix function $\frac{e^x - 1}{x}$ for non-invertible matrices I am trying to evaluate the matrix function
$$
f(X) = \frac{e^X - I}{X}, \tag{1}
$$
where $e^X$ is the usual matrix exponential, defined by
$$
e^X = \sum_{n=0}^{\infty} \frac{1}{n!} X^n.
$$
The problem is that all matrices I am working with are singular (if it is of relevance, $X \in \mathfrak{so}(3)$). According to this Wikipedia article, one may consider the function $f$ to be a real function, find its Maclaurin series, and use this to evaluate the function for a matrix. This seems reasonable enough, and it also seems to be how the matrix exponential is defined. The Maclaurin series for $(1)$ is given by
$$
f(x) = \sum_{n=0}^{\infty} \frac{X^n}{(n+1)!} = I + \frac{1}{2!}X + \frac{1}{3!}X^2 + \ldots,
$$
indicating that the invertibility of $X$ shouldn't matter when it comes to finding a meaningful value for $(1)$.
Now for my question:
Is there a closed form of $(1)$ that do not involve the inverse? To be clear, my goal is to implement the function in Python, and I would rather do it without computing a large partial sum and use that as the approximation. I know a closed form of the matrix exponential for matrices in $\mathfrak{so}(3)$.
 A: $\def\mymatrix#1{ \left[\begin{array}{c}#1\end{array}\right] }$
Nick Higham has a recommendation for calculating such matrix functions.
Write a simple function, applicable to non-singular matrices.
For example, here is some Julia code
function phi(X)
   return (exp(X)-I) / X
end

For every singular matrix there is a nearby non-singular matrix. So evaluate the simple function by perturbing the singular matrix with a random matrix $(X+R)\,$ where
$\,\|R\|\approx\sqrt{\varepsilon}\;$ where $\varepsilon$ is the machine epsilon. The function itself will be accurate to $\sqrt{\varepsilon}\approx10^{-8}$.
So for singular matrices, call the function like so
ε,X = 1e-8, rand(4,2)*rand(2,4)  # create a singular matrix
P = phi(X + ε*randn(4,4))        # evaluate the function
norm(I+X*P - exp(X))             # validate the function

If this simple method isn't good enough, there are rational (Pade) approximations for the exponential and phi-functions. Here is a fairly comprehensive paper.
Another suggestion is to evaluate the function using quad precision and a perturbation $\|R\|\approx10^{-16}.\;$ Then the function will be accurate to double precision (for singular matrix arguments).

Update
Another approach is to calculate the exponential of a block-triangular matrix
$$\eqalign{
\exp\left(\mymatrix{X&I\\0&0}\right)
 &= \mymatrix{e^X&\frac{e^X-e^0}{X}\\0&e^0}
 = \mymatrix{e^X&\frac{e^X-I}{X}\\0&I}
}$$
then extract the upper-right block. I'm not sure of the numerical accuracy of this method, but it should be as accurate as the underlying $\,\exp\,$ function. Note that this result holds for any analytic function, not just the exponential, i.e. 
$$\eqalign{
f\left(\mymatrix{X&I\\0&0}\right)
 &= \mymatrix{f(X)&\frac{f(X)-f(0)}{X}\\0&f(0)}
}$$
A: I suggest using the partial sum of power series expansion or
else a partial product of infinite product expansion. Here is
PARI/GP code for both methods:
/* pf(X) = pf(X/2)*(e^(X/2)+I)/2 */
pf(X, m=9) = {my(n=matsize(X)[1],
Id=matid(n), /* nxn identity matrix */
p=Id,        /* partial product matrix */
eX=Id+X/2^m  /* 1st order approx to e^(X/2^m) */); 
for(k=1,m, p*= (eX+1)/2; eX=sqr(eX)); p};

/* sf(X) = partial sum of power series */
sf(X, m=9) = sum(k=1,m, X^(k-1)/k!);

The partial sum of power series can be made more efficient by
using Horner's method. Try writing the equivalent code in Python
to see how the parameter m effects the accuracy of the result.
