Molecular vibrational frequency analysis plays an important role in theoretical and computational chemistry. However, in many cases, the analytical frequencies are unavailable, whereas frequency calculations using conventional numerical methods are very expensive. In this work, we propose an efficient method to numerically calculate the frequencies. Our main strategies are to exploit the sparseness of the Hessian matrix and to construct the N-fold two-variable potential energy surfaces to fit the parabola parameters, which are later used for the construction of Hessian matrices. A set of benchmark calculations is performed for typical molecules of different sizes and complexities using the proposed method. The obtained frequencies are compared to those calculated with the analytical methods and conventional numerical methods. It is shown that the results yielded with the new method are in very good agreement with corresponding accurate values (with a maximum error of ∼20 cm −1 ), while the required computation resource is largely reduced compared to that required by conventional numerical methods. For medium-sized molecules, the calculational scaling is lowered to O(N 1.6 ) (this work) from that of O(N 2 ) (conventional numerical methods). For even larger molecules, more computational savings can be achieved, and the scaling is estimated to be quasilinear with respect to the molecular size.