Optimizing a quadratic form integral over unit sphere

Question

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?

Answer 1

Since the $f_i(\theta,\phi)$ are linear combinations of spherical harmonics, we can write $$ \mathbf{f} = F \mathbf{Y} $$ where $\mathbf{Y}$ is a vector of the orthonormalized spherical harmonics - i.e.: $$ \int d\Omega Y_l^m Y_{l'}^{m'*} = \delta_{ll'} \delta_{mm'} $$ So the integral becomes, $$ \int d\Omega \mathbf{f}^{\dagger} (B^{\dagger} + B) \mathbf{f} = \int d\Omega \mathbf{Y}^{\dagger} F^{\dagger} (B^{\dagger} + B) F \mathbf{Y} $$ Note that since $\mathbf{Y}$ is orthonormal with respect to $\int d\Omega$, we have $$ \int d\Omega \mathbf{Y}^{\dagger} A \mathbf{Y} = Tr(A) $$ for any matrix $A$. This can be seen by explicitely writing out the integrand in components and then doing the trivial integrations using the orthonormality condition. So we have \begin{align} \int d\Omega \mathbf{f}^{\dagger} (B^{\dagger} + B) \mathbf{f} &= Tr(F^{\dagger} (B^{\dagger} + B) F) \\ &= 2 Tr(F^{\dagger} B F) \quad note: (Tr(A+A^{\dagger}) = 2 Tr(A)) \\ &= 2 Tr(Q L^{-1} F F^{\dagger}) \end{align} We want to find the $Q$ which maximizes this trace. By the answer given here, the solution is $$ Q = V U^{\dagger} $$ where $L^{-1} F F^{\dagger} = U S V^{\dagger}$ is the singular value decomposition of the matrix multiplying $Q$ in the trace.

Answer 2

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. 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.