The matrix function $\frac{e^x - 1}{x}$ for non-invertible matrices

Question

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)$.

Answer 1

$\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)} }$$

Answer 2

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.