The study of ground-state properties of the Fermi-Hubbard model is a long-lasting task in the research of strongly correlated systems, and is crucial for the understanding of notable quantum phenomena including superconductivity and magnetism. Owing to the exponentially growing complexity of the system, a quantitative analysis usually demands high computational cost and is restricted to small samples, especially in two or higher dimensions. Here, we introduce a variational method in the frame of fermionic Gaussian states, and obtain the ground states of one- and two-dimensional attractive Hubbard models via imaginary-time evolution. We calculate the total energy and benchmark the results in a wide range of interaction strength and filling factor with those obtained via exact two-body results, the density matrix renormalization group based on matrix product states (MPS), and projector Quantum Monte Carlo (QMC) method. For both 1D and 2D cases, the Gaussian variational method presents accurate results for total energy with a maximum systematic error ∼ 4% in the intermediate interaction region. The accuracy of these results has negligible dependence on the system size. We further calculate the double occupancy and find excellent agreement with MPS and QMC, as well as the experimental results of cold quantum gases in optical lattices. The results suggest that the Gaussian pairing state is a good approximation to the ground states of attractive Hubbard model, in particular in the strong and weak coupling limits. Moreover, we generalize the method to the attractive Hubbard model with a finite spin-polarization, which can be mapped to the repulsive interaction case via particle-hole transformation, and obtain accurate results of ground state energy and double occupancy. Our work demonstrates the ability of the Gaussian variational method to extract ground state properties of strongly correlated many-body systems with negligible computational cost, especially of large size and in higher dimensions.