# Optimizing a quadratic form integral over unit sphere

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?

• This is not a quadratic optimization problem. It is linear because $Q$ appears linearly in the objective function. It is only a well-posed problem because the set of orthogonal matrices is bounded. Aug 30, 2021 at 20:14
• You might find this answer from the math stack exchange helpful -- it looks like whoever posted that was trying to solve a similar problem. Aug 30, 2021 at 21:42
• is B dependent on the angular variables? Aug 31, 2021 at 15:02

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