Given two Hermitian matrices $\mathbf{A}$ and $\mathbf{B}$ of dimension $M \times M$, search for a set of $N \leq M$ orthonormal vectors $\mathbf{v}_n$ for which
$$\left( a_n \mathbf{A}+\mathbf{B} \right) \mathbf{v}_n=\lambda_n \mathbf{v}_n\;$$
where $\;a_n=\mathbf{v}_n^\dagger \mathbf{A} \mathbf{v}_n$, holds such that $\sum_n \lambda_n \rightarrow \max$.
Here, $\dagger$ denotes the complex conjugate transpose of a matrix/vector.
Is this problem solvable and if, how to solve it?
My thoughts so far:
First, the problem makes sense. Since $\mathbf{A}$ and $\mathbf{B}$ are hermitian both the $a_n$ and $\lambda_n$ are real and it makes sense to ask for the maximum of the sum.
My first idea was to solve the problem iteratively:
For each $n=1,\dots,N$ guess an initial $a_n$ (e.g. $a_n=1$), solve the eigenvalue problem and choose as $\mathbf{v}_n$ the eigenvector corresponding to the $n$-th largest eigenvalue $\lambda_n$, with this update $a_n$ and reiterate until self-consistency is reached. The problem with this method (besides possible convergence issues) is that each $\mathbf{v}_n$ comes from an individual diagonalization and thus the resulting $\mathbf{v}_n$ are not necessarily orthogonal.
The second idea/question that came to my mind was if the problem can be cast into a non-linear (quadratic?) eigenproblem that then can be solved with standard techniques?
Thank you very much for your suggestions, comments and advice!
P.S. (go on reading if you are interested in the origin of the problem)
This problem comes from an optimization problem in quantum physics. In the following, I will use the Bra-Ket notation commonly used in physics ($\langle f|g \rangle = \int \mathrm{d}\mathbf{r}\; f^\ast(\mathbf{r}) g(\mathbf{r})$ with $\ast$ denoting the complex conjugate). The goal is to minimize the functional $$\Omega[\lbrace w_n\rbrace ] = \sum_n [\langle w_n|\hat{\mathbf{R}}^2|w_n \rangle - \langle w_n|\hat{\mathbf{R}}|w_n \rangle^2]\rightarrow \mathrm{min}\;,$$ where $\hat{\mathbf{R}}$ is an operator and $\lbrace|w_n\rangle \rbrace$ a set of $N$ orthonormal states/functions ($\langle w_m|w_n \rangle = \delta_{mn}$). Taking the derivative $$\frac{\delta \Omega}{\delta w_n^\ast} = \hat{\mathbf{R}}^2 |w_n\rangle - 2\langle w_n|\hat{\mathbf{R}}|w_n \rangle \hat{\mathbf{R}} |w_n\rangle$$ and adding orthonormality as a constraint with Lagrangian multipliers $\Lambda_{mn}$, we need to solve $$\frac{\delta \Omega}{\delta w_n^\ast} + \sum_m \Lambda_{mn}\frac{\delta}{\delta w_n^\ast}[\langle w_n|w_m\rangle - \delta_{mn}] = 0\;.$$ Using a unitary transformation that diagonalizes $\Lambda$, this can be written as $$\frac{\delta \Omega}{\delta w_n^\ast} + \lambda_n |w_n\rangle = 0\;.$$ Now, expressing the states in an orthonormal basis $\lbrace|g_i\rangle\rbrace$ ($|w_n\rangle = \sum_{i=1}^M v_{in}|g_i\rangle$), this leads me to the above stated problem with the matrices $\mathbf{A}$ and $\mathbf{B}$ being the follwing matrix elements of the operators $A_{ij} = \sqrt{2}\langle g_i|\hat{\mathbf{R}}|g_j\rangle$ and $B_{ij} = -\langle g_i|\hat{\mathbf{R}}^2|g_j\rangle$.