Real symmetric matrix exponent

General matrix exponent is not trivial, scipy's implementation(help scipy.linalg.expm) refers a paper A New Scaling and Squaring Algorithm for the Matrix Exponential, and here is a revisited version.

But in most optimization cases, the matrices are real and symmetric, even positive definite. We can build a easy to understand method based on SVD. Let us assume that

$$A = USV'$$

is the SVD decomposition of a symmetric matrix A, we have

$$A = USV^T, A^T = VSU^T = A$$
$$A^2 = USV^TVSU^T = US^2U^T$$
$$A^3 = US^2U^TUSV^T = US^3V^T$$
$$A^{2k} = US^{2k}U^T$$
$$A^{2k+1} = US^{2k+1}V^T$$

This suggests that

$$e^A = \sum_k \dfrac{A^k}{k!} = \sum_k U\dfrac{S^{2k}}{(2k)!}U^T + \sum_k U\dfrac{S^{2k+1}}{(2k+1)!}V^T$$

Which leads to

$$e^A = \dfrac{U(e^S + e^{-S})U^T + U(e^S - e^{-S})V^T}{2}$$

If A is a positive definite matrix, then $U=V$, so the above formula can be simplified as $e^A=Ue^SU^T$.