Efficient way to solve generalized eigenvalue problem when the number of dimensions is greater than the number of samples

Question

I am trying to solve the generalized eigenvalue problem:

$$C_e v = \lambda C_o v.$$

$C_e$ and $C_o$ are both covariance matrices generated from data with $10512$ dimensions and about $2000$ samples. More specifically this is gridded forecast data, the grid has $10512$ points, and I have $2000$ samples. $C_e$ is generated from forecast data and $C_o$ is generated from analysis data. (This detail is inconsequential).

The dimensions of the covariance matrices, $C_e$ and $C_o$, are $10512\times 10512$.

The problem is solvable, but very computationally expensive for dimensions of those size.

Now, if this were just a standard eigenvalue problem, where I want to calculate the principal components of a data with $10512$ dimensions and $2000$ samples, the eigenproblem would be:

$$C_e v = \lambda v.$$

However, there is a trick where instead of using the covariance matrix, $C = 1/(N-1) \cdot X^\top X$, I could use the gram matrix, $G = 1/(N-1) \cdot X X^T$. The gram matrix will have dimensions of $2000 \times 2000$. Calculating the eigenvectors of this matrix is less computationally expensive than calculating the eigenvectors of $C$, which is $10512 \times 10512$. I could then use the eigenvectors of $G$ to calculate the eigenvectors of $C$. This is explained well here: Is PCA still done via the eigendecomposition of the covariance matrix when dimensionality is larger than the number of observations?

So, I am wondering if there is a similar "trick" to finding the eigenvectors for the generalize eigenvalue problem. I.e. maybe I can use the matrices $G_e$ and $G_o$ somehow?