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 .