I have an optimization problem, which is to maximize the following integral over the unit sphere: $$ \max_B \int d\Omega \mathbf{f}^{\dagger}(\theta,\phi) (B^{\dagger} + B) \mathbf{f}(\theta,\phi) $$ subject to: $$ B = Q L^{-1} $$ where $Q$ is an orthogonal matrix, and $L$ is a given lower triangular matrix. The vector $\mathbf{f} = (f_1, \dots, f_n) \in \mathbb{C}^n$ is known. So I need to find the orthogonal matrix $Q$ which maximizes the integral.
The functions $f_i(\theta,\phi)$ can be represented as a finite expansion in spherical harmonics, which means that given a $Q$, the integral can be evaluated exactly using Lebedev quadrature. The only method I can think of to solve this, is a Monte-Carlo approach where I generate random orthogonal matrices $Q$, calculate the integral with Lebedev quadrature, and save the $Q$ which maximizes the integral over some large number of $N$ trials. This seems quite inefficient, and may not yield a good $Q$ in the end anyway.
Can anyone think of a better approach to this optimization problem?