Here I basically do what joshphysics has already mentioned, just in a little more detail, and in a bit more intuitive basis (which makes effectively no difference). So definitely not a slick way.
I use the basis $|m_1\rangle\otimes|m_2\rangle$, where $m_i=\pm1/2$. Keeping only the sign we denote the basis as $\{|++\rangle,|+-\rangle,|-+\rangle,|--\rangle\}$.
Now for $S_z=S_{1z}+S_{2z}$ you can check that
$$
S_z|++\rangle=\left(\frac{\hbar}{2}+\frac{\hbar}{2}\right)|++\rangle
$$
and similarly
$$
S_z|--\rangle=\left(-\frac{\hbar}{2}-\frac{\hbar}{2}\right)|--\rangle
$$
so that you and up with the matrix
$$
S_z\rightarrow\hbar\left[\begin{array}{cccc}1&0&0&0\\0&0&0&0\\0&0&0&0\\0&0&0&-1\end{array}\right]
$$
and
$$
S_z^2\rightarrow\hbar^2\left[\begin{array}{cccc}1&0&0&0\\0&0&0&0\\0&0&0&0\\0&0&0&1\end{array}\right]
$$
For $S_x$ and $S_y$ you use raising and lowering operators as josh suggested and obtain the matrix elements by acting on each state in the basis and seeing what has a non-zero dot product with the result, e.g.
$$
\begin{align}
S_y|+-\rangle=\frac{1}{2i}\left(S_+-S_-\right)|+-\rangle&=\\\frac{1}{2i}\left(S_{1+}+S_{2+}-S_{1-}-S_{2-}\right)|+-\rangle&=\frac{1}{2i}\left(\hbar|++\rangle-\hbar|--\rangle\right)
\end{align}
$$
thus $\langle++|S_y|+-\rangle=\frac{\hbar}{2i}$ and $\langle--|S_y|+-\rangle=-\frac{\hbar}{2i}$. Once you do this you can get
$$
S_x\rightarrow\hbar\left[\begin{array}{cccc}0&1&1&0\\1&0&0&1\\1&0&0&1\\0&1&1&0\end{array}\right]
$$
and
$$
S_y\rightarrow\hbar\left[\begin{array}{cccc}0&1&1&0\\-1&0&0&1\\-1&0&0&1\\0&-1&-1&0\end{array}\right]
$$
Finally the Hamiltonian looks like
$$
H=\left[\begin{array}{cccc}\hbar^2B&0&0&\hbar^2A\\0&0&0&0\\0&0&0&0\\\hbar^2A&0&0&\hbar^2B\end{array}\right]
$$
Therefore $|+-\rangle$ and $|-+\rangle$ are eigenvectors with eigenvalue $0$, and you only need diagonalize $2\times2$ matrix to find that the eigenvalue $\hbar^2\left(B+A\right)$ corresponds to the eigenvector $\left(|++\rangle+|--\rangle\right)/\sqrt{2}$, and $\hbar^2\left(B-A\right)$ corresponds to the eigenvector $\left(|++\rangle-|--\rangle\right)/\sqrt{2}$.
The set $G$ gives the representation of the identity and generators of the abstract group of quaternions as elements in $SL(2,\mathbb C)$ which are also in $SU(2)$. Taking the completion of this yields the representation $Q_8$ of the quaternions presented in the question.
From the description of the symmetry group as coming from here, consider the composition of two $\pi$ rotations along the $\hat x$, $\hat y$, or $\hat z$ axis. This operation is not the identity operation on spins (that requires a $4\pi$ rotation). However, all elements of $D_2$ given above are of order 2.
This indicates that the symmetry group of the system should be isomorphic to the quaternions and $Q_8$ is the appropriate representation acting on spin states. The notation arising there for $D_2$ is probably from the dicyclic group of order $4\times 2=8$ which is isomorphic to the quaternions.
Best Answer
The commutator of an observable like $S_x$ or $P_y$ (the spin in the x direction and the y component of the momentum respectively) with the hamiltonian will tell you about the time evolution of that observable. It tells you how that observable changes with time.
This is given through the heisenberg equation:
$\frac{dA}{dt} = \frac{i}{\hbar}[H, A]$
where A is an operator (like $S_x$ or $P_y$) and I have assumed that A does not depend explicitly on time, i.e $A = A(x(t),p(t))$ and not $A = A(x(t), p(t), t)$.
Now from the equation above, we can see that if A commutes with the hamiltonian, $[H,A] = 0$, then $\frac{dA}{dt} = 0$ and thus A is constant in time, it is conserved.
So when you showed that $S_z$ commutes with H while $S_x$ and $S_y$ do not, you showed that the z component of a particles spin is constant in time, while the x and y components are not constant, they precess.