As shown in the image, I have a 16QAM constellation that is misaligned due to a phase rotation. In this case you can see that rotation amount is approximately $\theta = \pi/4$, but this won't be the case in general. For real-world data the phase could be a slowly varying function of time, $\theta(t)$, so that it is not enough to apply some fixed correction factor.
I am aware of differential mapping schemes that solve the phase ambiguity problem due to the constellation having $\pi/2$ symmetry, but it seems $\theta$ must still be known to perform the slicing.
One suggested solution was to try to map the received constellation point to the nearest QAM constellation point and feed a phase-locked-loop with the result, but it is not clear how this would perform when $\theta$ varies over time.
What techniques exist to recover the symbols? I have already tried various carrier recovery schemes based on feedback loops, with no success, and am interested in decision directed approaches that may avoid having to find the phase.