Traditional methods for characterizing an optimized molecular structure as a minimum or as a saddle point on the nuclear potential energy surface require the full Hessian. However, if f denotes the number of nuclear degrees of freedom, a full Hessian calculation is more expensive than a single point geometry optimization step by the order of magnitude of f. Here we present a method which allows to determine the lowest vibrational frequencies of a molecule at significantly lower cost. Our approach takes advantage of the fact that only a few perturbed first-order wave functions need to be computed in an iterative diagonalization scheme instead of f ones in a full Hessian calculation. We outline an implementation for Hartree-Fock and density functional methods. Applications indicate a scaling similar to that of a single point energy or gradient calculation, but with a larger prefactor. Depending on the number of soft vibrational modes, the iterative method becomes effective for systems with more than 30-50 atoms.