Gaussian Process Factor Analysis (GPFA) hasbeen broadly applied to the problem of identi-fying smooth, low-dimensional temporal struc-ture underlying large-scale neural recordings.However, spike trains are non-Gaussian, whichmotivates combining GPFA with discrete ob-servation models for binned spike count data.The drawback to this approach is that GPFApriors are not conjugate to count model like-lihoods, which makes inference challenging.Here we address this obstacle by introduc-ing a fast, approximate inference method fornon-conjugate GPFA models. Our approachuses orthogonal second-order polynomials toapproximate the nonlinear terms in the non-conjugate log-likelihood, resulting in a methodwe refer to aspolynomial approximate log-likelihood(PAL) estimators. This approxima-tion allows for accurate closed-form evalua-tion of marginal likelihoods and fast numericaloptimization for parameters and hyperparam-eters. We derive PAL estimators for GPFAmodels with binomial, Poisson, and negativebinomial observations. We find the PAL esti-mation achieves faster convergence times andhigh accuracy compared to existing state-of-the-art inference methods. We also find thatPAL hyperparameters can provide sensible ini-tialization for black box variational inference(BBVI), which will improve BBVI accuracy.We apply these methods to data from mousevisual cortex and primate higher-order visualand parietal cortices. We demonstrate thatPreliminary work. Under review by AISTATS 2020. Do notdistribute.PAL estimators achieve fast and accurate ex-traction of latent structure from multi-neuronspike train data.