We explore the mantle density structure of the northeast Atlantic region using constrained linear inversion of the satellite gravity gradient data based on statistical prior information and assuming a Gaussian model. The uncertainty of residual gravity gradient signal is characterized by covariance matrix obtained using geostatistical analysis of controlled-source seismic data. The forward modelling of the gravity gradients in the 3D reference crustal model is performed using a global spherical harmonics analysis. We estimate the model covariance function in the radial and angular directions using a variogram method. We compute volumetric gravity gradient kernels for a spherical shell covering the northeast Atlantic region down to the upper limit of the mantle transition zone (410 km depth). The solution of the linear inverse problem in the form the mean density model follows a least-squares approach. The results indicate a direct relationship between the seismic velocity and density anomalies in the Iceland-Jan Mayen region, Greenland and the Norwegian passive margin. The predicted low-density anomalies at the depth of 100-150 km underneath the northeast Atlantic Ocean are correlated with the distribution of Cenozoic underwater volcanoes and seamount-like features of the seafloor.