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 M1,2 and calculating the corresponding tidal deformabilities Λ1,2, we perform Monte Carlo integration over the GW likelihood for M1,2 and Λ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 Λ1.4 = 160 +448 −113 and p(2ρnuc) = 1.35 +1.8 −1.2 × 10 34 dyn/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 Λ1.4 = 556 +163 −172 and p(2ρnuc) = 4.73 +1.4 −2.5 × 10 34 dyn/cm 2 . We further infer the maximum NS mass supported by the EOS to be Mmax = 2.09 +0.37 −0.16(2.04 +0.22 −0.002 ) M with the loose (tight) prior. The Bayes factor between the two priors is B A I 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. *