In the context of seismic Full Waveform Inversion(FWI), the Hessian matrix is used to determine the model resolution. Having an exact computation of the Hessian matrix will lead to more reliable inferences involving parameters of the subsurface. One of the main problems however is its high computational cost, making the use of the Hessian matrix not practical, particularly in large offset experiments. In this work, we propose the estimation of the Hessian matrix from a sub-sampled set of its columns. The proposed approach computes a subset of Hessian-vector products using the Second Order Adjoint-State Method(SOASM), as proposed by Metivier et al. (2013), and then estimates the Hessian matrix using the compressed sampling(CS)theory. Our approach computes only a set of linear combinations of the columns of the matrix and then estimates the whole Hessian matrix by using the Gradient Projection for Sparse Reconstruction(GPSR)algorithm. The reconstruction of the Hessian matrix is feasible since it can be considered a sparse matrix. This novel approach leads to better results than approximation methods such as the Broyden-Fletcher-Goldfarb-Shannon(BFGS)algorithm, and it requires less computation time than computing the full matrix with the SOASM method. We tested the proposed approach by reconstructing a Hessian matrix from a number of Hessian-vector products equivalent to 60% of the columns of the full matrix. The reconstructed Hessian matrix using the CS method resulted in a PSNR of 53.3 dB when it is compared with the original Hessian matrix. In contrast, the Hessian matrix obtained with the BFGS method results in a PSNR of 28.1 dB.