How to prove that any rotation can be represented by 3 Euler angles How can one prove that any rotation of a rigid object in 3-dimensional (3D) space can be represented by a sequence of three rotations around pre-fixed axes by 3 Euler angles? I see this statement in many textbooks, but so far I did not find a proof of the statement. 
I do understand that 3 parameters are generally needed to represent rotation of a 3D object (e.g. from Goldstein, Poole, and Safko, Classical Mechanics, 3rd Ed. Ch. 4). However, I can not be sure that 3 Euler angles can be such 3 parameters.  
For now, I accept that, for any 3D rotation, there exists a unique matrix in the group SO(3) that transforms the coordinates of a point in the rotated object. I also accept that such a matrix for the rotation around the $z$-axis by angle $\alpha$ is expressed as
\begin{equation}
R_z (\alpha) = 
\begin{pmatrix}
\cos \alpha  & \sin \alpha & 0 \\
-\sin \alpha & \cos \alpha & 0 \\
0 & 0 & 1
\end{pmatrix},
\end{equation}
and that rotation around $y$-axis has a similar matrix representation.
Therefore, I understand that for a given $A \in \mathrm{SO}(3)$, I can calculate the Euler angles by comparing 3 elements of $A$ (e.g. some of the elements in the 3rd row and 3rd columns in the convention adopted below) with the following product of 3 matrices corresponding to the rotations by 3 Euler angles,
\begin{equation}
R_z (\alpha)
R_y (\beta)
R_z (\gamma) = 
\begin{pmatrix}
\cos \alpha  & \sin \alpha & 0 \\
-\sin \alpha & \cos \alpha & 0 \\
0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
\cos \beta  & 0 & -\sin \beta \\
0 & 1 & 0 \\
\sin \beta & 0 & \cos\beta
\end{pmatrix}
\begin{pmatrix}
\cos \gamma  & \sin \gamma & 0 \\
-\sin \gamma & \cos \gamma & 0 \\
0 & 0 & 1
\end{pmatrix}.
\end{equation}
However, I can not be sure if the calculated Euler angles always equate the other elements not used in the calculation.
Is there a simple proof without going through much algebraic manipulation?
 A: An algorithm that solves for $\alpha$, $\beta$, and $\gamma$ for any given proper 3×3 rotation matrix constitutes a constructive proof.

\begin{equation}
R_z (\alpha)
R_y (\beta)
R_z (\gamma) = 
\begin{pmatrix}
\cos \alpha  & \sin \alpha & 0 \\
-\sin \alpha & \cos \alpha & 0 \\
0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
\cos \beta  & 0 & -\sin \beta \\
0 & 1 & 0 \\
\sin \beta & 0 & \cos\beta
\end{pmatrix}
\begin{pmatrix}
\cos \gamma  & \sin \gamma & 0 \\
-\sin \gamma & \cos \gamma & 0 \\
0 & 0 & 1
\end{pmatrix}
\end{equation}

Multiplying this out yields
\begin{equation}
R_{zyz} = R_z (\alpha)
R_y (\beta)
R_z (\gamma) =
\begin{pmatrix}
\text{hot mess} & \text{hot mess} & -\cos\alpha \sin\beta \\
\text{hot mess} & \text{hot mess} & \phantom{+}\sin\alpha \sin\beta \\
\sin\beta \cos\gamma & \sin\beta \sin\gamma & \cos\beta
\end{pmatrix}
\end{equation}
Note that the last element of the last row is $\cos\beta$. Given some proper rotation matrix $R$, one can find $\beta$ by taking the inverse cosine of that element:
\begin{equation}
\beta = \arccos(R_{3,3})
\end{equation}
Note that the other two elements of the last column yield an expression for $\alpha$, and that the other two elements of the last row yield an expression for $\gamma$:
\begin{equation}
\begin{aligned}
\alpha &= \arctan\left(\frac{\phantom{+}R_{2,3}}{-R_{1,3}}\right) \\
\beta &= \arctan\left(\frac{\phantom{+}R_{3,2}}{\phantom{+}R_{3,1}}\right)
\end{aligned}
\end{equation}
You'll need to use the two argument version of $\arctan$ to disambiguate between (for example) 60 degrees (pi/3) and 240 degrees (4/3*pi). I'm not going to go through the math, but this solution replicates the "hot mess" that I ignored.
There's one issue with this solution, and that's when $\sin\beta = 0$. This happens when $R_{3,3}$ is ±1. The remaining elements in the last column and last row are identically zero in this case. This is called "gimbal lock". My "hot mess" becomes a lot simpler in this case:
\begin{equation}
R_{zyz} =
\begin{pmatrix}
\pm\cos(\alpha+\gamma)  & \pm\sin(\alpha+\gamma) & 0 \\
\mp\sin(\alpha+\gamma) & \mp\cos(\alpha+\gamma) & 0 \\
0 & 0 & \pm 1
\end{pmatrix}
\end{equation}
You can solve for $\alpha+\gamma$ in the case of gimbal lock, but there is no unique solution for $\alpha$ and $\gamma$. That's okay; just choose an arbitrary value (typically zero) for $\gamma$ (or $\alpha$) and solve for the other given that you can solve for $\alpha+\gamma$.

Your zyz sequence is rather odd, in two regards. The canonical Euler rotation sequence is a rotation about z followed by a second rotation about the once-rotated x axis followed third rotation about the twice-rotated z axis. Yours is a rotation about the initial z axis, followed by a second rotation about the initial y axis, followed by a third rotation about the initial z axis. Those are but two of the twenty four different rotation sequences that are oftentimes called Euler angles. In all twenty four cases, you'll find that


*

*Exactly one element of the rotation matrix is either $\sin \beta$, $-\sin\beta$, or $\cos \beta$, where $\beta$ represents the second rotation in the sequence (which lets you solve for $\beta$),

*That the other four elements on the same row / same column as this special element are not a "hot mess",

*Except in the case of gimbal lock, these yields unique solutions for the other two elements of the rotation sequence, and

*In the case of gimbal lock, you can solve for $\alpha+\gamma$ or $\alpha-\gamma$. In that case, arbitrarily choose a value for $\gamma$ (typically zero) and you'll be able to find $\alpha$.

A: There's a constructive proof that can be understood intuitively. I assume z is vertical and y is forward/backward. You rotate the object about the z-axis until the top is somewhere on the xz plane i.e. y=0. This makes the top of the object point perpendicular to the y-axis, so you can rotate about the y-axis until it points up. And now you just need to rotate it on the z-axis until forward points in the right direction.
If you want to know what the original rotation was, just take the opposite of each of those steps and put them in reverse order. For example, if you rotated 10°, -30°, 50°, then to get the original rotation from the correct rotation it's just -50°, 30°, -10°.
This sequence of operations can always be done regardless of how an object is oriented. It's not necessarily unique. If up is already pointing up, then it will be on the correct plane no matter how much you rotate about the z-axis. But the point is that there is some rotation about the z-axis that leaves it on the correct plane, which there is.
Edit:
If you want something strictly more mathematical, suppose you have an orthogonal frame of vectors, $x = (x_1,x_2,x_3), y = (y_1,y_2,y_3),$ and $z = (z_1,z_2,z_3)$.
$z$ is the top of the object, so first we rotate it so $z'_2 = 0$. We just rotate it by $\arctan\frac{z_2}{z_1}$, and it ends up pointing at $z' = \left(\pm\sqrt{z_1^2+z_2^2},0,z_3\right)$. And rotate it another $180^\circ$ if that got a negative instead of a positive, so it's $z' = \left(\sqrt{z_1^2+z_2^2},0,z_3\right)$. $\arctan\frac{z_2}{z_1}$ is only really undefined in the case of $z_1=z_2=0$, in which case don't rotate it at all.
Now we rotate on the $y$-axis by $\arctan\frac{\sqrt{z_1^2+z_2^2}}{z_3}$ and we get $z'' = \left(0,0,\sqrt{z_1^2+z_2^2+z_3^2}\right)$ which must be $(0,0,1)$ since it's a unit vector. Again, if it's undefined, we don't need to rotate it.
Since $y \perp z$ and we're only rotating, $y'' \perp z''$.
From there, we know $y''_3 = 0$, so we have $y'' = (y''_1,y''_2,0)$.
Just rotate it on the $z$-axis by $\arctan\frac{y''_1}{y''_2}$ plus an extra $180^\circ$ if it faces the wrong way, and we get $y''' = \left(0,\sqrt{y''^2_1+y''^2_2},0\right)$. And that's just $(0,1,0)$, since it's a unit vector.
In this case $y''_1$ and $y''_2$ can't both be zero, so we don't have to worry about $\arctan$ being undefined at all.
Since we rotated on the $z$-axis, and $z''$ was already on the $z$-axis, $z''' = z''$.
All we have left is $x$. Since $x = y \times z$, and we're only rotating, $x''' = y''' \times z'''$.
$x''' = (0,1,0) \times (0,0,1)$
$x''' = (1,0,0)$
A: For convenience, I would first like to change the sign convention. That is,
\begin{equation}
R_{z}(\alpha) = \begin{pmatrix}\cos\alpha&-\sin\alpha&0\\\sin\alpha&\cos\alpha&0\\0&0&1\end{pmatrix},\quad
R_{y}(\beta) = \begin{pmatrix}\cos\beta&0&\sin\beta\\0&1&0\\-\sin\beta&0&\cos\beta\end{pmatrix},\quad
R_{z}(\gamma) = \begin{pmatrix}\cos\gamma&-\sin\gamma&0\\\sin\gamma&\cos\gamma&0\\0&0&1\end{pmatrix},
\end{equation}
and we consider
\begin{equation}
R(\alpha,\beta,\gamma) \equiv R_{z}(\alpha) R_{y}(\beta) R_{z}(\gamma).
\end{equation}
Next, recall that the 3D rotation matrix $R(\alpha,\beta,\gamma)$ is uniquely determined by specifying where the three unit vectors $\hat{e}_{1}\equiv(1,0,0)^{T}$, $\hat{e}_{2}\equiv(0,1,0)^{T}$, and $\hat{e}_{3}\equiv(0,0,1)^{T}$ are mapped to. [These are in fact the three columns of $R(\alpha,\beta,\gamma)$.]
It is straightforward to check that $\hat{e}_{3}$ is mapped to
\begin{equation}
\hat{e}_{3}^{\prime}\equiv R(\alpha,\beta,\gamma)\,\hat{e}_{3} =(\sin\beta\cos\alpha,\sin\beta\sin\alpha,\cos\beta)^T
\equiv \hat{n}(\beta,\alpha),
\end{equation}
where $\hat{n}(\beta,\alpha)$ is the unit vector with the polar angle $\beta$ and the azimuthal angle $\alpha$. Observe that by tweaking $\alpha$ and $\beta$, we can make $R(\alpha,\beta,\gamma)\,\hat{e}_{3}$ to be any unit vector we choose.
Next, suppose that we fix $\alpha$ and $\beta$, and set $\gamma=0$. We have $\hat{e}_{3}^{\prime} = \hat{n}(\beta,\alpha)$, and the two unit vectors $\hat{e}_{1}^{\prime}$ and $\hat{e}_{2}^{\prime}$ (defined analogously to $\hat{e}_{3}^{\prime}$) satisfy the following conditions:
(1) Both lie on the plane that is perpendicular to $\hat{n}(\beta,\alpha)$ and contains the origin.
(2) They have a fixed orientation [because $\mathrm{det}R(\alpha,\beta,\gamma)=1$].
(3) $e_{1}^{\prime} \perp e_{2}^{\prime}$.
By further rotating about $\hat{n}(\beta,\alpha)$ by an arbitrary angle $\gamma$, we can adjust $\hat{e}_{1}^{\prime}$ and $\hat{e}_{2}^{\prime}$ to be any two unit vectors satisfying the constraints stated above. Hence, an arbitrary 3D rotation matrix $M$ can be represented by
\begin{equation}
M = R_{\hat{n}(\alpha,\beta)}(\gamma)R(\alpha,\beta,0) = R_{\hat{n}(\alpha,\beta)}(\gamma) R_{z}(\alpha)R_{y}(\beta).
\end{equation}
Here, $\alpha$ and $\beta$ determine $M \hat{e}_{3}$, and for a fixed $M \hat{e}_{3}$, $\gamma$ determines $M \hat{e}_{1}$ and $M \hat{e}_{2}$.
Now notice that $R(\alpha,\beta,0)$ maps the $z$ unit vector (i.e., $\hat{e}_{3}$) to $\hat{n}(\alpha,\beta)$. Hence,
\begin{equation}
 R_{\hat{n}(\alpha,\beta)}(\gamma) = R(\alpha,\beta,0) R_{z}(\gamma) R(\alpha,\beta,0)^{-1},
\end{equation}
and
\begin{equation}
M = R_{\hat{n}(\alpha,\beta)}(\gamma)R(\alpha,\beta,0) = R(\alpha,\beta,0) R_{z}(\gamma) = R_{z}(\alpha)R_{y}(\beta) R_{z}(\gamma) = R(\alpha,\beta,\gamma).
\end{equation}
