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](https://en.wikipedia.org/wiki/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?