Approximating invariant subspaces of generalized eigenvalue problems (GEPs) is a fundamental computational problem at the core of machine learning and scientific computing. It is, for example, the root of Principal Component Analysis (PCA) for dimensionality reduction, data visualization, and noise filtering, and of Density Functional Theory (DFT), arguably the most popular method to calculate the electronic structure of materials. Given Hermitian $H,S\in\mathbb{C}^{n\times n}$, where $S$ is positive-definite, let $Π_k$ be the true spectral projector on the invariant subspace that is associated with the $k$ smallest (or largest) eigenvalues of the GEP $HC=SCΛ$, for some $k\in[n]$. We show that we can compute a matrix $\widetildeΠ_k$ such that $\lVertΠ_k-\widetildeΠ_k\rVert_2\leq ε$, in $O\left( n^{ω+η}\mathrm{polylog}(n,ε^{-1},κ(S),\mathrm{gap}_k^{-1}) \right)$ bit operations in the floating point model, for some $ε\in(0,1)$, with probability $1-1/n$. Here, $η>0$ is arbitrarily small, $ω\lesssim 2.372$ is the matrix multiplication exponent, $κ(S)=\lVert S\rVert_2\lVert S^{-1}\rVert_2$, and $\mathrm{gap}_k$ is the gap between eigenvalues $k$ and $k+1$. To achieve such provable "forward-error" guarantees, our methods rely on a new $O(n^{ω+η})$ stability analysis for the Cholesky factorization, and a smoothed analysis for computing spectral gaps, which can be of independent interest. Ultimately, we obtain new matrix multiplication-type bit complexity upper bounds for PCA problems, including classical PCA and (randomized) low-rank approximation.