What is the link between the rotating wave approximation and the algebraic representation of a dynamical system?

I think it may be easier to understand the system of equations if you swap the ordering of operators in ($\star$) to $(a,b,a^*,b^*)$. In this case the equation becomes

$$ \frac{d}{dt} \left( \begin{array}{c} a \\ b \\ a^* \\ b^* \end{array} \right) = -i H_{\textrm{eff}} \left( \begin{array}{c} a \\ b\\ a^* \\ b^* \end{array} \right) \tag{$\star$} \, . $$ $$ \begin{align*} H_{\textrm{eff}} & = \sigma_z \otimes \left( g \sigma_x + \frac{d}{2} \sigma_z + \frac{s}{2} \mathbb{I} \right) - i\, g (\sigma_y \otimes \sigma_x) \\ &= \left( \begin{array}{cccc} \omega_1 & g & 0 & -g \\ g & \omega_2 & -g & 0 \\ 0 & g & -\omega_1 & -g \\ g & 0 & -g & -\omega_2 \end{array} \right) \end{align*} $$ The $ i\, g\, \sigma_y\otimes \sigma_x $ term then corresponds to the block off-diagonal part of the matrix. The first diagonal block has eigenvalues $\frac{s}{2} \pm \frac{1}{2}\sqrt{d^2 + 4 g^2}$, where $s = \omega_1 + \omega_2, d = \omega_1-\omega_2$. To get the second block's eigenvalues we send $s\rightarrow-s$. So for $s\gg |d|, g$, these blocks are split in 'energy' by approximately $s$.

Once the block off-diagonal terms are removed, each of the two remaining blocks contains only both positive or both negative frequencies, i.e. the clockwise or counterclockwise rotating parts of each physical mode. Therefore, one intuitive way to understand the meaning of the rotating wave approximation is that it keeps only dynamical terms that rotate in the same direction.

The 'secular approximation' for this 'Hamiltonian' corresponds to throwing away terms that are 1) small and 2) do not 'conserve energy'. In the case $s\gg |d|, g$, the term $ i\, g\, \sigma_y\otimes \sigma_x $ is small and changes the energy of a 'state' by approximately $\pm s$. This is unlike the other term proportional to $g$, which although is small does not change the energy of a 'state' by more than $\sim |g|$.