The virtual two-loop corrections for Higgs production in gluon fusion are calculated analytically in QCD for arbitrary Higgs and quark masses. Both scalar and pseudo-scalar Higgs bosons are considered. The results are obtained by expanding the known one-dimensional integral representation in terms of mH/mq, and matching it with a suitably chosen ansatz of Harmonic Polylogarithms. This ansatz is motivated by the known analytic result for the Higgs decay rate into two photons. The method also allows us to check this result and to extend it to the pseudo-scalar decay rate.