I am seeking a clever solution to the following problem. Given $$X(\theta,\phi) = exp(-iA(\theta,\phi))\; B\; exp(+iA(\theta,\phi))$$ with the square, Hermitian matrix $A$: $$A(\theta,\phi) = A_{0,0} Y_{0,0}(\theta,\phi) + A_{1,0} Y_{1,0}(\theta,\phi) + A_{1,1} Y_{1,1}(\theta,\phi) + A_{1,-1} Y_{1,-1}(\theta,\phi)$$ and Hermitian $X(\theta,\phi)$, I want to determine the decomposition of $X(\theta,\phi)$ in spherical harmonics, i.e., $$X(\theta,\phi) = \sum_{l,m} X_{l,m} Y(\theta,\phi).$$ Currently, I am evaluating $X(\theta,\phi)$ on a grid of angles and determine the expansion coefficients by the method of least squares based on the discretized $X(\theta_i,\phi_i)$.
Is there a more elegant approach that would lead me to $X_{l,m}$ directly (in particular for large $l$)?
The $A_{l,m}$s and $B$ do not commute, $B$ is a Hermitian matrix independent of $\theta$ and $\phi$, $i$ is the imaginary unit, and the $Y_{l,m}$s are spherical harmonics.