Efficient Algorithms for the Matrix Cosine and Sine

Several improvements are made to an algorithm of Higham and Smith for computing the matrix cosine. The original algorithm scales the matrix by a power of 2 to bring the $\infty$-norm to 1 or less, evaluates the [8/8] Pad\'e approximant, then uses the double-angle formula $\cos(2A) = 2\cos^2A - I$ to recover the cosine of the original matrix. The first improvement is to phrase truncation error bounds in terms of $\norm{A^2}^{1/2}$ instead of the (no smaller and potentially much larger quantity) $\norm{A}$. The second is to choose the degree of the Pad\'e approximant to minimize the computational cost subject to achieving a desired truncation error. A third improvement is to use an absolute, rather than relative, error criterion in the choice of Pad\'e approximant; this allows the use of higher degree approximants without worsening an a priori error bound. Our theory and experiments show that each of these modifications brings a reduction in computational cost. Moreover, because the modifications tend to reduce the number of double-angle steps they usually result in a more accurate computed cosine in floating point arithmetic. We also derive an algorithm for computing both $\cos(A)$ and $\sin(A)$, by adapting the ideas developed for the cosine and intertwining the cosine and sine double angle recurrences