Assuming $B$ is not angle-dependent, you know the coefficients of each $f_i$ a priori, and $\mathbf{f}$ can only be evaluated by summing up the spherical harmonics expansion, simplify your problem by making use of the orthogonality of spherical harmonics instead of numerically integrating.
Let $f_i(\theta,\phi)=\sum_{lm}f_{i,lm}Y_l^m(\theta,\phi)$. Then your cost function becomes, using the summation convention: $$ \begin{align} g &= \int d\Omega f_{i,lm}^Y_l^{m}(\theta,\phi) (B_{ji}^ + B_{ij})f_{j,l'm'}Y_{l'}^{m'}(\theta,\phi)\ &= f_{i,lm}^f_{j,l'm'}(B_{ji}^ + B_{ij}) \int d\Omega Y_l^{m}(\theta,\phi) Y_{l'}^{m'}(\theta,\phi)\ &= f_{i,lm}^f_{j,l'm'}(B_{ji}^ + B_{ij}) \delta_{ll',mm'}\ &= f_{i,lm}^f_{j,lm}(B_{ji}^ + B_{ij}) \end{align} $$ Assuming $F_{lm}=(f_{1,lm},\ldots,f_{n,lm})$, this can be written as $$ g=\sum_{lm} F_{lm}^\dagger (B^\dagger+B) F_{lm}. $$ We see that for each scalar element of the sum $F_{lm}^\dagger B^\dagger F_{lm}=(F_{lm}^\dagger B F_{lm})^\dagger$, so we can simplify further to $$ g=2\mathfrak{R}\sum_{lm} F_{lm}^\dagger Q L^{-1} F_{lm}. $$ Since you only have control over the elements of $Q$, you can precompute and store as a vector $x_{lm}=L^{-1}F_{lm}$s. 804
Now, you know the $F_{lm}$ and $x_{lm}$ coefficients, and $Q$ is by definition orthogonal. One way to make the sum large is to make $Q$ span the space of the "largest" $F_{lm},x_{lm}$. You could take the $k \leq n$ largest-norm $F_{lm}$ and $x_{lm}$ vectors, put them in a matrix (horzcat in Matlab), take the SVD, and use its column space matrix as your $Q$.
A similar approach can be easily derived in the case that $\mathbf{f}$ is cheap to evaluate (i.e., with fewer than $\mathcal{O}(L^2)$ operations per quadrature point, where $L$ is the maximal degree of spherical harmonics used), with samples of $\Omega$ taking the place of the $l,m$ indices. I realize that taking the real part in the above expression may change your approach slightly; have not yet thought that through completely. 976