Algorithms for computing matrix logarithm. On my quest to find the holy grail of mathematics become a little bit better at algebra, I have read up on matrix logarithms and exponentials and how useful they can be in investigating groups and algebras. I know that power methods are common for trancendental matrix functions. This question is about fast numerical ways to compute the matrix logarithm.

Own work: I know of and have implemented the Taylor expansion: $$\log({\bf I + A}) = \sum_{i=1}^N\frac{(-1)^i}{i}{\bf A}^i$$ or equivalently $$\log({\bf A}) = \sum_{i=1}^N\frac{(-1)^i}{i}({\bf A-I})^i$$ in various ways. Having stored $({\bf A-I})^{i-1}$ at iteration $i$ lets us multiply with $\bf (A-I)$ just once every iteration. I have read about some popular speed-up technique using the fact that $$2\log({\bf A}^{1/2}) = \log({\bf A})$$ and that Taylor series is more accurate close to the point of expansion for the lower exponents. However as I don't know any fast matrix-square root, I have not yet managed to employ this fact. Maybe a Taylor expansion of the square root function would do?
 A: Power methods are actually not so common (anymore?) for computing matrix functions, especially if your matrix has a high condition number.
Try Padé approximants, inverse scaling-and-squaring, or a Schur algorithm instead.  The matrix logarithm and the matrix square root are well explored topics in the field of matrix function computations, so it's best to just check out the established sources.


*

*Chapter 11 of Functions of Matrices by Nicholas Higham.  Nick Higham is THE expert on matrix functions, so pretty much anything he's written about them is gold.

*If you can't get a hold of Higham's books, try one of his survey articles.  Section 5.2 pertains to the matrix logarithm in this article.

*I haven't read this one, but it appears to be another good source.

A: 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}$.
