Computing the logarithm of a symmetric positive definite matrix (Q1294459)

The following algorithm to compute \(\log(A)\) for a matrix \(A>0\) is proposed. (1) Reduce to tridiagonal form: \(A=QTQ^T\), (2) Compute an approximant \[ R_m(X)=\sum_{j=1}^m a_j (I+b_jX)^{-1}X \] with \(X=\mu T-I\), (3) Set \(S_m=-\log\mu I+QR_mQ^T\). The approximant \(R_m(x)\) is a diagonal Padé approximant for \(\log(1+x)\) with \(x\in(-1,1)\), so that the \(a_j\) and \(b_j\) are the weights and abscissas of the \(m\)-point Gauss-Legendre quadrature formula. The degree \(m\) is selected such that a certain precision is obtained. This \(m\) (and also the parameter \(\mu\)) can be computed in terms of the largest and smallest eigenvalue of \(A\). An easy estimate for the optimal \(m\) in function of the condition number of \(A\) is also derived from the error estimate of the Padé approximant. The complexity of the algorithm is analysed and several illustrative numerical examples are included.