  Since their first discovery, quasars have been essential probes of the
distant Universe. However, due to our limited knowledge of its nature,
predicting the intrinsic quasar continua has bottlenecked their usage. Existing
methods of quasar continuum recovery often rely on a limited number of
high-quality quasar spectra, which might not capture the full diversity of the
quasar population. In this study, we propose an unsupervised probabilistic
model, Quasar Factor Analysis (QFA), which combines factor analysis (FA) with
physical priors of the intergalactic medium (IGM) to overcome these
limitations. QFA captures the posterior distribution of quasar continua through
generatively modeling quasar spectra. We demonstrate that QFA can achieve the
state-of-the-art performance, $\sim 2\%$ relative error, for continuum
prediction in the Ly$\alpha$ forest region compared to previous methods. We
further fit 90,678 $2<\mathrm{z}<3.5$, SNR$>2$ quasar spectra from Sloan
Digital Sky Survey Data Release 16 and found that for $\sim 30\%$ quasar
spectra where the continua were ill-determined with previous methods, QFA
yields visually more plausible continua. QFA also attains $\lesssim 1\%$ error
in the 1D Ly$\alpha$ power spectrum measurements at $\mathrm{z}\sim 3$ and
$\sim 4\%$ in $\mathrm{z}\sim 2.4$. In addition, QFA determines latent factors
representing more physically motivated than PCA. We investigate the evolution
of the latent factors and report no significant redshift or luminosity
dependency except for the Baldwin effect. The generative nature of QFA also
enables outlier detection robustly; we showed that QFA is effective in
selecting outlying quasar spectra, including damped Ly$\alpha$ systems and
potential Type II quasar spectra.
