Algorithms for computing matrix logarithm.

Let $A\in M_n(\mathbb{C})$ be a matrix that has no eigenvalues in $\mathbb{R}^{-}$. We want to calculate an approximation of $\log(A)$, the principal log of $A$.

There is (eventually) an instability (in the calculation) when the eigenvalues of $A$ are close or when the matrix of change of basis (to obtain a diagonal form, otherwise a Jordan form) has a large condition number. Note that, when we use a numerical method, we cannot know if $A$ is diagonalizable or not (especially if we know only an approximation of $A$). However, it is often the case that eigenvalues ​​are very close. You can use the algorithms exposed by Higham but they are usually complicated. When $A$ is not too mean, you can try that

Step 1. You calculate approximations of the eigenvalues of $A$ (these approximations are always distinct) and you calculate $p$, the Lagrange interpolation polynomial between $(\lambda_i)_i$ and $(\log(\lambda_i))_i$. Then you calculate $B=p(A)$ and you have an approximation of $L=\log(A)$ (the complexity of calculation is $O(n^4))$.

Step 2. $B-I_n+Ae^{-B}$ is a better approximation.

Indeed $e^B-A=e^B-e^L\approx (B-L)e^B$ because $L,B$ commute (they are both polynomials in $A$). Thus $L-B\approx -I_n+Ae^{-B}$. Note that $e^{-B}$ is easy to calculate (divide $B$ by $2^k$).

Example: $A$ is similar to $diag(2,4,2.00001,1.99999,3.99999)$; then $||A-e^B||\approx 10^{-6}$ and $||A-\exp(B-I+Ae^{-B})||\approx 4.10^{-14}$.

One simple way is as follows. The principal logarithm of a matrix whose eigenvalues do not lie on the closed negative real axis can be written as \begin{equation*} \log\boldsymbol{A}=\int_0^1 \left[(1-\eta)\boldsymbol{I}+\eta\boldsymbol{A}\right]^{-1}(\boldsymbol{A}-\boldsymbol{I})\,d\eta. \end{equation*} One can convert the above integral to standard Gaussian quadrature by using the substitution \begin{equation*} \eta=\frac{1+\xi}{2}, \end{equation*} so that we have \begin{equation*} \log\boldsymbol{A}=\int_{-1}^1\left[(1-\xi)\boldsymbol{I}+(1+\xi)\boldsymbol{A}\right]^{-1}(\boldsymbol{A}-\boldsymbol{I})\,d\xi. \end{equation*} Now use standard Gaussian quadrature to compute the above integral (you will find programs on the net that generate weights and locations of the Gauss points for any order of quadrature). You can download a fortran implementation matrixlog.f and data file matrixlog.dat from https://mecheng.iisc.ac.in/csjog/Codes/