Glacial isostatic adjustment (GIA) is the response of the solid Earth to past ice loading, primarily, since the Last Glacial Maximum, about 20 K yrs BP. Modeling GIA is challenging because of large uncertainties in ice loading history and also the viscosity of the upper and lower mantle. GPS data contain the signature of GIA in their uplift rates but these also contain other sources of vertical land motion (VLM) such as tectonics, human and natural influences on water storage that can mask the underlying GIA signal. In this study, we use about 4000 GPS vertical velocities as observational estimates of global GIA uplift rates, after correcting for major elastic deformation effects. A novel fully-automatic strategy is developed to post-process the GPS time series and to correct for non-GIA artefacts. Before estimating vertical velocities and uncertainties, we detect outliers and jumps and correct for atmospheric mass loading displacements. We correct the resulting velocities for the elastic response of the solid Earth to global changes in ice sheets, glaciers, and ocean loading, as well as for changes in the Earth's rotational pole relative to the 20th century average. We then apply a spatial median filter to remove sites where local effects are dominant to leave approximately 4000 GPS sites. The resulting novel global GPS dataset shows a clean GIA signal at all post-processed stations and is therefore suitable to investigate the behavior of global GIA forward models. The results are transformed from a frame with its origin in the center of mass of the total Earth's system (CM) into a frame with its origin in the center of mass of the solid Earth (CE) before comparison with 13 global GIA forward model solutions, with best fits with Pur-6-VM5 and ICE-6G predictions. The largest discrepancies for all models were identified for Antarctica and Greenland, which may be due to either uncertain mantle rheology, ice loading history/magnitude and/or GPS errors.