Modeling quasar spectra is a fundamental task in astrophysics as quasars are the tell-tale sign of cosmic evolution. We introduce a novel unsupervised learning algorithm, Quasar Factor Analysis (QFA), for recovering the intrinsic quasar continua from noisy quasar spectra. QFA assumes that the Lyα forest can be approximated as a Gaussian process, and the continuum can be well described as a latent factor model. We show that QFA can learn, through unsupervised learning and directly from the quasar spectra, the quasar continua and Lyα forest simultaneously. Compared to previous methods, QFA achieves state-of-the art performance for quasar continuum prediction robustly but without the need for predefined training continua. In addition, the generative and probabilistic nature of QFA paves the way to understanding the evolution of black holes as well as performing out-of-distribution detection and other Bayesian downstream inferences.