I've heard some people mistakenly claim that the Moon does not rotate around an axis (further claiming that this is the reason we always see the same hemisphere; evidently it's not), but here is a purely mathematical argument that it almost certainly does rotate around an axis (i.e. with probability one) using only linear algebra:
Step 1: First let's make the reasonable simplifying assumptions that the Moon is a perfect $3$-dimensional ball $B^3$ of unit radius ($1$ unit $=$ radius of the moon), and that the space it's moving in is Euclidean 3-space.
Step 2: Next let's make sure to distinguish between tranlational motion in which the ball's center actually changes its position in space, and stationary motion where the position of the center point remains fixed while the rest of the points on the ball may change position--the latter is the type of motion we are interested in. It seems to me at this point before going through any argument and with all things being equally likely, the probability that the Moon is undergoing some stationary motion is equal to one (there's only one way to be still, but there are infinitely many ways to move). It only remains to decide what that stationary motion is.
Step 3: Next let's argue that we may use linear transformations to model the sphere's stationary motion: Neglecting any translational motion, if we imagine taking snapshots of our ball at time $t_0$ and again at time $t_1>t_0$, then the transformation of the $t_0$ points of our ball to the $t_1$ points defines a function
$$F\colon B^3\rightarrow B^3.$$
Further this function should satisfy the following properties:
It fixes the center point, which we may identify as the origin, $0$, i.e.
$$F(0)=0.$$
It fixes distances which we may express in terms of inner products on $\mathbb{R}^3$ as
$$(F(v),F(w))=(v,w), \ \forall \ v,w\in B^3.$$
If fixes orientations, which we may express in terms of determinants on $\mathbb{R}^3$ as
$$F(\mathbf{e_1})\wedge F(\mathbf{e_2})\wedge F(\mathbf{e_3})>0$$
where $\mathbf{e_1},\mathbf{e_2},\mathbf{e_3}$ is the standard basis in $\mathbb{R}^3$.
To justify items 2. and 3., imagine fixing three rigid rods (of some fixed lengths, all less than $1$ unit) to the center of the Moon pointing in linearly independent directions--say the standard basis directions--at time $t_0$. Then at time $t_1$ the rods may be pointing in different directions, but they should still have the same lengths, and they should still satisfy the same right hand rule, i.e. have same orientation.
First note that we can always extend $F\colon B^3\rightarrow B^3$ to a function
$$\hat{F}\colon\mathbb{R}^3\rightarrow\mathbb{R}^3$$
by
$$\hat{F}(v)=\begin{cases} F(v) & \text{if $v\in B^3$}\\ |v|\cdot F\left(\frac{v}{|v|}\right) & \text{if $v\notin B^3$}\\ \end{cases}.$$
One can further check that $\hat{F}$ also satisfies conditions 1, 2, and 3 above.
Exercise: Show that $\hat{F}\colon {\mathbb{R}}^3\rightarrow\mathbb{R}^3$ is actually a linear transformation.
Hint: Compute the inner products
$$(\hat{F}(v_1+v_2)-\hat{F}(v_1)-\hat{F}(v_2),\hat{F}(v_1+v_2)-\hat{F}(v_1)-\hat{F}(v_2))$$
$$(\hat{F}(cv)-c\hat{F}(v),\hat{F}(cv)-c\hat{F}(v))$$
and use property 2.
Step 4: Now that we have a linear function $\hat{F}\colon\mathbb{R}^3\rightarrow\mathbb{R}^3$ modeling the stationary motion of the Moon from time $t_0$ to $t_1$, we can compute its standard matrix $A$. This will be a $3\times 3$ real orthogonal matrix (from Property 2) with $\det(A)=1$ (by Property 3). Its characteristic polynomial $c_A(\lambda)$ will be a cubic polynomial with real coefficients. Since such polynomials always have at least one real root, we find that $A$ must have at least one real eigenvalue.
Exercise: Show that one of these real eigenvalues must equal $1$.
(Of course the eigenspace for $\lambda=1$ will turn out to be exactly the axis of rotation for our Moon)
Step 5: If $\ell\in \mathbb{R}^3$ is the $A$-eigenspace for $\lambda=1$, let $W\subset\mathbb{R}^3$ be its orthogonal compliment. Then of course $\hat{F}$ restricts to give a linear transformation
$$\hat{F}'\colon W\rightarrow W.$$
Choosing any basis $\mathcal{B}'$ for $W$, we can compute the $\mathcal{B}'$-matrix for $\hat{F}'$, say $A'$. Then of course $A'$ will be a real $2\times 2$ orthogonal matrix with $\det(A')=1$.
Exercise: Show that every real $2\times 2$ orthogonal matrix $A'$ with $\det(A')=1$ is a rotation through some angle $\theta$, i.e. is similar to a matrix of the form
$$A'=\left(\begin{array}{rr} \cos(\theta) & -\sin(\theta)\\ \sin(\theta) & \cos(\theta)\\ \end{array}\right).$$
Now of course this means that we can represent our function $\hat{F}\colon\mathbb{R}^3\rightarrow\mathbb{R}^3$ by a $3\times 3$ matrix in block diagonal form as
$$A=\left(\begin{array}{rrr} 1 & 0 & 0\\ 0 & \cos(\theta) & -\sin(\theta)\\ 0 & \sin(\theta) & \cos(\theta)\\ \end{array}\right).$$
In particular the stationary motion taking $B^3$ at time $t_0$ into $B^3$ at time $t_1$ was rotation through an angle $\theta$ around the axis $\ell$.
Step 6: I suppose that if we make the additional simplifying assumption that the stationary motion of the moon is also continuous as a function of time, then fixing $t_0$ and letting $t_1=t$ vary over the time interval $[t_0,\infty)$ we get a continuous $1$-parameter family of linear transformations
$$t\mapsto\left\{\hat{F}_t\colon \mathbb{R}^3\rightarrow\mathbb{R}^3\right\}.$$
Then since the roots of a (characteristic) polynomial are continuous in its coefficients, we could deduce that the axis $\ell=\ell_t$ (i.e. eigenspace for $\lambda(t)\equiv 1$) remains fixed for all time and hence so does its perp space $W=W_t$, so our continuous motion can really be modeled by the $1$-parameter family of matrices
$$t\mapsto A_t=\left(\begin{array}{rrr} 1 & 0 & 0\\ 0 & \cos(\theta(t)) & -\sin(\theta(t))\\ 0 & \sin(\theta(t)) & \cos(\theta(t))\\ \end{array}\right)$$
presumably where the $\theta(t)$ is a continuous function in $t$.
I presented this argument (which I think is due to the French mathematician Chasles) to my linear algebra students last year, and I think they found it interesting :-)
