Understanding the operation of Thomas precession How can we physically understand the operation of Thomas Precession? 
Thomas Precession is an effect that modifies the effective energy of coupling between the spin and the orbital angular momentum of the electron by an extra factor of $\frac{1}{2}$. I know that Thomas precession is related to the fact that two boosts in different directions do not commute and is equivalent to a pure boost plus a pure rotation. 
But how does this lead to a precessional motion of the spin?
How does this effect physically modifies the dynamics of the electron inside the atom such that the interaction energy picks up an extra factor of $\frac{1}{2}$?
 A: It's a little unclear what you are asking, but not to worry: I find the hardest part of working something out is finding the right question in the first place. So I'll answer two questions I perceive in the above and we can go from there.
You seem to be asking:


*

*For an account of Thomas precession wholly in terms of the Wigner rotation (the rotation part of the composition of two noncommuting boosts);

*For further insight into the Wigner rotation itself;
For the first question, I think one is meant to imagine the electron as a classical point particle and work out out what happens to its spin as it "orbits" the nucleus with relativistic effects accounted form. Then you "guess" the Hamiltonian from the relativistically corrected classical interaction. Take heed that the missing factor of $1/2$ is not a problem with the Dirac-equation-described electron, which is why I've never given this problem much thought. 
Okay, so we have three inertial frames of reference:


*

*That of the nucleus;

*That which is momentarily co-moving with the electron in a circular orbit, which is say velocity $v$ in the $x$-direction relative to the first frame, when the electron path's tangent is in this direction;

*That of the electron a short time $\delta t$ later. Relative to the second frame, this frame will be moving at velocity $a_y \hat{y} \delta t + a_x \hat{x} \delta t$, where I have split a general acceleration into components parallel and at right angles to the velocity of frame 2 relative to frame 1 (although $a_x=0$ for circular motion).


To get the Lorentz transformation mapping from frame 1 to 3 directly, we compose boosts:
$$\exp\left(\frac{a_y}{c} \,\delta t\, K_y + \frac{a_x}{c} \,\delta t\, K_x\right) \exp\left(\frac{v_x}{c}\,K_x\right) = \exp\left(\frac{v_x}{c}\,K_x + \frac{a_x}{c}\,\delta t\,K_x + \frac{a_y}{c}\,\delta t\, K_y + \frac{a_y\,v_x}{2\,c^2}\,[K_y,\,K_x]\,\delta t + \cdots\right)\tag{1}$$ 
by the Dynkin formula version of the Baker-Campbell-Hausdorff theorem. From the commutation relationships for the Lorentz group we get $[K_y,\,K_x] = i\,J_z$, so that $\frac{a_y\,v_x}{2\,c^2}\,[K_y,\,K_x]\,\delta t$ corresponds to a rotation about the $z$ axis of $a_y\,v_x\,\delta t / (2\,c^2)$ radians (here, somewhat obviously $J_z,\,J_y,\,J_z$ are "generators of rotations", i.e. tangents to rotations about the $x,\,y,\,z$ axes at the identity, respectively, whilst $J_z,\,J_y,\,J_z$ are "generators" of boosts).
So, for circular motion, $a_x = 0$ and we find, from the Dynkin formula that the rate of spin precession is:
$$\vec{\Omega} =  \frac{\vec{a}\wedge \vec{v}}{2\,c^2} + \cdots\tag{2}$$
whence we see the factor of $1/2$ which then enters the dynamical equations used to get the classical interaction expression Equation 29 in the reference[1] below that gives the right quantum Hamiltonian when used to "guess" the latter. Actually, there are other, smaller terms in the CBH formula also of order $\delta t$ that sum up to:
$$\vec{\Omega} =  \frac{\gamma^2\,\vec{a}\wedge \vec{v}}{(\gamma+1)\,c^2}\approx \frac{\vec{a}\wedge \vec{v}}{2\,c^2}\tag{3}$$
(the approximation holding for $v\ll c$) but the Dynkin formula lets you see the largest contribution to the Wigner rotation very clearly, and gives the "right" factor which yields the "right" classical interaction. Reference [2] below derives the full Wigner rotation expression from the Rodrigues formula for $SL(2,\,\mathbb{C})$, which is the double cover of the group of proper, orthochronous Lorentz transformations (the identity-connected component of the Lorentz group). Also see referece [3], which derives the commutator discussed above by an elementary trick, but it is still only accurate to first order in $a / c$. 
[1]: G. F. Smoot,UC Berkeley Physics 139 Lecture notes, 1998;
[2]: Dr. Bill Pezzaglia, CSUEB Physics, "Special Relativity
and Thomas Precession", 2010, page 16 onwards.
[3]: Andrzej Dragan, Tomasz Odrzygóźdź, "Half-page derivation of the Thomas precession", 2012, arxiv.org/abs/1211.1854
Now, for my second "perceived question", if you want further insight into the Wigner rotation itself, then I don't believe one can get a deeper understanding than your statement:

...two boosts in different directions do not commute and is equivalent to a pure boost plus a pure rotation

The language of Lie groups is the "right way" to describe these ideas insofar it is the most concisely accurate one we have so far. There is no "everyday" insight. It simply arises from the structure constants of the Lorentz algebra, and for me this is as deep as it gets. Although boosts in one direction are a one-parameter subgroup of the Lorentz group, the set of all boosts is not a subgroup of the Lorentz group, and this fact is simply owing to the nature of the Lorentz group. Now, there may well be, even probably will be, other descriptions of the above in the future, but I'd bet they'll be thought up by a mathematician as an alternative to the Lorentz algebra commutation relationships. The effect cannot be understood in everyday terms, it is wholly relativistic and indeed even the ratio:
$$\frac{\frac{a_y\,v_x}{2\,c^2}\,[K_y,\,K_x]\,\delta t}{\frac{a_x}{c}\,\delta t\,K_x + \frac{a_y}{c}\,\delta t\, K_y}\tag{4}$$
of largest non-Galilean term in the Dynkin formula $\frac{a_y\,v_x}{2\,c^2}\,[K_y,\,K_x]\,\delta t$ to the Galilean change in velocity $\frac{a_x}{c}\,\delta t\,K_x + \frac{a_y}{c}\,\delta t\, K_y$ vanishes in the Galilean limit ($c\to\infty$). 
To get a feel for the depth of your statement, try describing the following analogous, everyday phenomenon precisely in words other than the language you have used in your statement, i.e. the languages of commutators and Lie groups. The following is a phenomenon and question that is utterly everyday and that will be well wonted to anyone who lives in Mumbai, Tokyo, Beijing, Paris, Los Angeles or Sydney and who has driven a car: 
How do we parallel-park a car, i.e. induce a pure translation orthogonal to the car's heading, when we can only go forwards and backwards on curved paths?
I don't believe you can get a better (more everyday, concise) description that is also precise than:
Steer operators of different path curvatures do not commute and indeed always compose to a steer operator that is composed with a pure translation.
If we represent our car's configuration as a $2\times1$ column vector of two complex numbers $(z,\,u)^T$ where $z$ stands for the position of the car and $u=\exp(i\,\theta)$ encodes the car's orientation. Then, with an Ackermann steering mechanism, if we drive the car a distance $s$ along a path of constant curvature (the steering wheel set to a constant setting), we shall evolve the car's configuration following:
$$\mathrm{d}_s\left(\begin{array}{c}z\\u\end{array}\right)=\left(\begin{array}{cc}0&1\\0&i\,\kappa\end{array}\right)\left(\begin{array}{c}z\\u\end{array}\right)\tag{5}$$
where $\kappa\in\mathbb{R}$ encodes the steering setting as the curvature of the path; it can be any setting in some nonzero interval $-\kappa_0\leq\kappa\leq +\kappa_0$ where $\kappa_0^{-1}$ is radius of the car's tightest turning circle. The Lie algebra of the set of all possible transformations on the car is defined by:
$$\begin{array}{lclcl}
D &=& \text{"drive"}&=& \left(\begin{array}{cc}0&1\\0&0\end{array}\right)\\
R &=& \text{"rotate"}&=& \left(\begin{array}{cc}0&0\\0&i\end{array}\right)\\
S &=& \text{"sideslide"}&=& \left(\begin{array}{cc}0&i\\0&0\end{array}\right)\\
\end{array}\tag{6}$$
with the commutation relationships:
$$[D,\,R]=S;\quad [S,\,D] = 0;\quad [R,\,S] = D\tag{7}$$
so that the problem is wholly analogous to the Thomas precession idea: all the transformations available to us are of the form $\exp((D + \kappa R) s)$ where $s$ is the distance driven, and such transformations always compose to give a third transformation of the same form, composed with a pure translation (the analogue here of the Wigner rotation), which is what we need (a pure, orthogonal translation) to get out of the parallel carpark. This ends the answer, but ....

Background on the Carparking Analogue of the Wigner Rotation
To complete my discussion of this wonderful little problem, if we compose three manoeuvres of the kind $\exp((D + \kappa R) s)$ that the car's steering affords us, i.e. we impart the transformation $\exp((D+\kappa_3 R)s_3)\exp((D+\kappa_2 R)s_2)\exp((D + \kappa_1 R)s_1)$, and if we also ensure the condition $s_1\kappa_1 + s_2\kappa_2+s_3\kappa_3 = 0$, we get the pure translation represented by the addition of the complex number:
$$i\,(\kappa_1-\kappa_2) \left(1-e^{i\,s_1\,\kappa_1}\right)+i\,(\kappa_2-\kappa_3) \left(1-e^{i\,s_3\,\kappa_3}\right)\tag{8}$$
to the car's $z$ co-ordinate. Moreover, this can be purely imaginary (i.e. a pure "sideslide") for arbitrarily small $s_1,\,s_2,\,s_3$ if we ensure:
$$(\kappa_1-\kappa_2) \sin(s_1\,\kappa_1)+(\kappa_2-\kappa_3) \sin(s_3\,\kappa_3) = 0\tag{9}$$
Since $s_1,\,s_2,\,s_3$ can be arbitrarily small, we can do this in an arbitrarily tight park. However, (8) means that the translation is second order in the $s_j$, so that the number of times we must repeat this sequence of three manoeuvres is proportional to $w/(L\,\epsilon^2)$ , where the parking space’s length is $1+\epsilon$ times the length $L$ of the car, whose width is $w$. Indeed, with our basic manoeuvres of the form $\exp((D + \kappa R) s)$ we can realise any transformation in the group:
$$\begin{array}{l}\begin{array}{lcl}E_2 &=& \left\{\mathcal{E}(x,\,y,\,\theta):\; x,\,y,\,\theta\in\mathbb{R}\right\}\\
\mathcal{E}(x,\,y,\,\theta)&\stackrel{def}{=}&
\exp\left(\left(\begin{array}{cc}
 0 & (i\,x-y) \mathcal{A}(\theta) \\
 0 & i \theta\end{array}\right)\right) = \left(
\begin{array}{cc}
 1 & x+i\,y \\
 0 & e^{i \theta}\\
\end{array}
\right)\\
\mathcal{A}(\theta) &\stackrel{def}{=}& \left\{\begin{array}{lll} \frac{\theta}{e^{i\,\theta} - 1}&;&\theta\neq 0\\1&;&\theta = 0\end{array}\right.\\
\end{array}\\
\mathcal{E}(x_2,\,y_2,\,\theta_2) \cdot\mathcal{E}(x_1,\,y_1,\,\theta_1)= \mathcal{E}(x_1 + e^{i\,\theta_1} x_2,\,y_1 + e^{i\,\theta_1} y_2,\,\theta_1+\theta_2)\end{array}\tag{10}$$
which is the whole Euclidean group of affine transformations of the Argand plane. Howeveer, even though we can realise any member of this group as a nett transformation, we can't realise the continuous path of the form $\{e^{s\,R}:\,s\in\mathbb{R}\}$ or of the form $\{e^{s\,S}:\,s\in\mathbb{R}\}$; we can only reach any given point on this path by way of a zigzaggy path keeping arbitrarily near to these continuous paths.
