  We develop a non-parametric method for inferring the universal neutron star
(NS) equation of state (EOS) from gravitational wave (GW) observations. Many
different possible realizations of the EOS are generated with a Gaussian
process conditioned on a set of nuclear-theoretic models. These synthetic EOSs
are causal and thermodynamically stable by construction, span a broad region of
the pressure-density plane, and can be selected to satisfy astrophysical
constraints on the NS mass. Associating every synthetic EOS with a pair of
component masses $M_{1,2}$ and calculating the corresponding tidal
deformabilities $\Lambda_{1,2}$, we perform Monte Carlo integration over the GW
likelihood for $M_{1,2}$ and $\Lambda_{1,2}$ to directly infer a posterior
process for the NS EOS. We first demonstrate that the method can accurately
recover an injected GW signal, and subsequently use it to analyze data from
GW170817, finding a canonical deformability of $\Lambda_{1.4} =
160^{+448}_{-113}$ and $p(2\rho_{\mathrm{nuc}})=1.35^{+1.8}_{-1.2}\times
10^{34}~\mathrm{dyn}/\mathrm{cm}^2$ for the pressure at twice the nuclear
saturation density at 90$\%$ confidence, in agreement with previous studies,
when assuming a loose EOS prior. With a prior more tightly constrained to
resemble the theoretical EOS models, we recover $\Lambda_{1.4} =
556^{+163}_{-172}$ and $p(2\rho_{\mathrm{nuc}})=4.73^{+1.4}_{-2.5}\times
10^{34}~\mathrm{dyn}/\mathrm{cm}^2$. We further infer the maximum NS mass
supported by the EOS to be $M_\mathrm{max}=2.09^{+0.37}_{-0.16}$
($2.04^{+0.22}_{-0.002}$) $M_\odot$ with the loose (tight) prior. The Bayes
factor between the two priors is $B^{\mathcal{A}}_{\mathcal{I}} \simeq 1.12$,
implying that neither is strongly preferred by the data and suggesting that
constraints on the EOS from GW170817 alone may be relatively prior-dominated.
