Can the exponential map of $SO(3)$ be expressed as a limit? According to this textbook (page 14), the exponential map of some Lie algebras (like $SO(2)$) can be thought of as repeated multiplication of matrices with infinitesimally small rotations.
$$
\lim_{k\to\infty} \left(\mathbf{I} + \alpha_a\frac{X_a}{k}\right)^k = \sum_{m=0}^{\infty} \frac{1}{m!}\left(\alpha_aX_a\right)^m \equiv e^{\alpha_aX_a},
$$
where $X_a$ is the $a$-th generator (basis vector) of the Lie algebra and $\alpha_a$ is the coefficient of the $a$-th basis vector (the textbook seems to be using Einstein's summation notation).
I feel like this derivation of the exponential map for "sphere-like" manifolds like $SO(2)$ or $SO(3)$ makes sense and thus I wanted to test whether this intuitive definition via a limit, when simplified, can indeed lead to the analytical version of the exponential map.
I therefore first decided to use Python's symbolic math package SymPy to verify whether the equation holds for the Lie algebra of $SO(2)$ and SymPy was able to simplify the symbolic limit into the analytical form of the exponential map (3.6) in the textbook. I plugged the generator $X_{SO(2)}$ into the limit and after simplification, I got the expected analytical form $\exp_{SO(2)}(x)$.
$$
\begin{aligned}
X_{SO(2)} &= \begin{bmatrix}0 & -x\\x & 0\end{bmatrix}\\
\exp_{SO(2)}(x) &= \begin{bmatrix}\cos x & -\sin x\\\sin x & \cos x\end{bmatrix}
\end{aligned} 
$$
However, when I tried to simplify the limit with generators for $SO(3)$, the analytical expression and the limit were not the same!
I wondered whether the same limit could be used to derive the analytical form of the exponential map of $SO(3)$, so I found the generator $X_{SO(3)}$ of $SO(3)$ and the expected analytical exponential map $\exp_{SO(3)}(X)$ in this paper (page 5, Example 3 and 4):
$$
\begin{aligned}
X_{SO(3)} &= \begin{bmatrix}0 & -z & y\\ z & 0 & -x\\ -y & x & 0\end{bmatrix}\\
\exp_{SO(3)}(X) &= \mathbf{I} + X_{SO(3)} \sin{\theta} + X^2_{SO(3)} (1 - \cos{\theta})
\end{aligned} 
$$
I assume that $\theta = \sqrt{x^2+y^2+z^2}$.
I expected the following to be true:
$$
\lim_{k\to\infty} \left(\mathbf{I} + \frac{X_{SO(3)}}{k}\right)^k \stackrel{?}{=} \mathbf{I} + X_{SO(3)} \sin{\theta} + X^2_{SO(3)} (1 - \cos{\theta})\\
\lim_{k\to\infty} \left(\mathbf{I} + \frac{X_{SO(3)}}{k}\right)^k - \left(\mathbf{I} + X_{SO(3)} \sin{\theta} + X^2_{SO(3)} (1 - \cos{\theta})\right) \stackrel{?}{=} \begin{bmatrix}0 & 0 &0\\0 & 0 &0\\0 & 0 &0\end{bmatrix}
$$
However, after attempting to use SymPy to simplify the expression above, I did not get a $\mathbf{0}$ matrix. It seems like one of the following is true:

*

*I might be using a wrong analytical form $\exp_{SO(3)}(x,y,z)$.

*I might misunderstood what $\theta$ stands for in $\exp_{SO(3)}(x,y,z)$.

*I made a mistake when re-writing the expressions into SymPy.

*SymPy is unable to correctly simplify the expression.

*The limit does not hold for $SO(3)$.

What did I do wrong? And if the limit does not hold for $SO(3)$, what are some other ways to derive the exponential map for $SO(3)$ in an intuitive and satisfying way?
Below is the Python/SymPy code I used:
from sympy import *
init_printing()


x, y, z = symbols('x y z')
G = Matrix([
    [0, -z, y],
    [z, 0, -x],
    [-y, x, 0]
])

# The analytic expression of the exponential map
omega = sqrt(x*x + y*y + z*z)
analytic_expression = Identity(3).as_explicit() + sin(omega) * G + (1 - cos(omega)) * G**2

# Expressing the exponential map as a limit
n = symbols('n', integer=True)
limit_expression = (Identity(3).as_explicit() + G/n)**n
limit_expression = limit_expression.limit(n, oo)

# If both expressions are exactly the same, their difference should be the zero matrix
diff_matrix = simplify(analytic_expression - limit_expression)

# diff_matrix is NOT zero. Why?
print(diff_matrix)


EDIT:
To rule-out the possibility of a bug in SymPy, I have also tried to simplify the expression using Wolfram Mathematica (I assume Mathematica is more reliable than SymPy). However, not even Mathematica simplified the difference between the limit and the analytical solution into the $\mathbf{0}$ matrix. The Mathematica code:
Clear["Global`*"]
X = {
    {0, -z, y},
    {z, 0, -x},
    {-y, x, 0}
}

omega = Sqrt[x^2 + y^2 + z^2]

analyticExpression = IdentityMatrix[3] + Sin[omega]*X + (1 - Cos[omega]) * (X.X)
limitExpression = Limit[MatrixPower[IdentityMatrix[3] + X/n, n], n -> Infinity]

(analyticExpression - limitExpression) // FullSimplify // MatrixForm

 A: The problem is $e^{X}=I + Y \sin\theta + Y^2 (1-\cos\theta)$ where $Y=\frac{1}{\theta} X$, not $X$ itself.
The identity can be proved by diagonalizing $X$:
$$\lim_{k\rightarrow\infty} (I+\frac{X}{k})^k = I+Y\sin\theta + Y^2(1-\cos\theta)$$
Note that the matrix $X$ really just represents the linear map $\vec v\mapsto (x, y, z)^t\times \vec v$. It's obvious that $(x,y,z)^t$ is an eigenvector of eigenvalue $1$ for both sides. And both sides rotates the orthogonal complement of $(x,y,z)^t$ by angle $\theta$ (in other words, $X$ has $i\theta$ as an eigenvalue, and $Y$ has $i$ as an eigenvalue, and $e^X$ has $e^{i\theta}$ as an eigenvalue). These can be easily confirmed by the diagonalizing the matrix $X$.
Or to put this in a more self-contained and rigorous standing, it's easy to see that for $X\not=0$, it has three distinct eigenvalues $0, i\theta, -i\theta$ with corresponding eigenvectors $v_0, v_1, v_2$ separately. For each eigenvalue $\lambda$ and corresponding eigenvector $v$, $$\lim_{k\rightarrow\infty} (I+\frac{X}{k})^k v=\lim_{k\rightarrow \infty} (I+\frac{X}{k})^{k-1}(1+\frac{\lambda}{k})v = \cdots = \lim_{k\rightarrow\infty} (1+\frac{\lambda}{k})^kv=e^\lambda v$$
And it's easy to check that $$(I+Y\sin\theta + Y^2(1-\cos\theta))v_0=v_0$$ $$(I+Y\sin\theta + Y^2(1-\cos\theta))v_1=(\cos\theta + i\sin\theta)v_1=e^{i\theta}v_1$$ etc. Therefore, the two sides match when acting on a basis of $\mathbb C^3$, hence coincide as matrices.
Or you can just use the ODE theory to justify the (more general) result as done in the references.
A: Response to comment:
As the WP article provided, and your source!, insist, the generators have to be normalized!
$$
X_{SO(3)}= \theta K, ~~~~\hbox{s.t.} ~~ \operatorname{Tr}K^2=-2, \\
K^3=-K. $$
It is then self-explanatory (no computer algebra nonsense!) that
$$
e^X= e^{\theta K}= I+\sin\theta ~ K + (1-\cos\theta)~ K^2,
$$
the celebrated Rodrigues rotation formula, (1840).
Response to further comment
The point of k=3 was to help you appreciate the combinatorics and frame the general case in terms of the binomial theorem,
$$
\left ( I + {\theta K\over k}\right ) ^k = \sum_{n=0}^k  \frac{k!}{n!(k-n)! k^n} ~\theta^n K^n 
\\ = I + \left (\theta - \frac{\theta^3}{3!}(1+O(1/k))+...\right )K + \left ({\theta^2\over 2} (1+O(1/k)+...\right ) K^2,
$$
to appreciate your limit cannot fail! (mindful of $K^3=-K$,  $K^4=-K^2$, $K^5=K$, ... )
