Influence functional theory has been applied to describe vibrational energy relaxation of molecules in the solution based upon harmonic oscillators bath approximation. An algebraic formula of the perturbative influence functional is presented for a number of quantum bath oscillators (q(k)) nonlinearly coupled to the system x as Sigma(k)f(k)(x)q(k) + Sigma(k)Sigma(l)g(kl)(x)q(k)q(l) + Sigma(k)Sigma(l)Sigma(m)h(klm)(x)q(k)q(l)q(m). The approach opens a way to molecular based analysis of multiphonon processes making usage of a number of techniques and concepts in the field of path integral and quantum field theory. Based upon the functional, we also derive a computationally tractable expression for the relaxation time by executing the path integral exactly. The theory is of much higher approximation than Fermi＇s golden rule including perturbations up to the infinite order. A recipe for the numerical work based upon classical molecular dynamics calculation followed by the normal mode analysis is also presented.