The matrix function $\frac{e^x - 1}{x}$ for non-invertible matrices
$\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).