We present a complete analytical expression of the two-dimensional (2D) vibrational response function without invoking the factorization approximations based on the Wick's theorem that have been used in most previously reported theories. Since the harmonic approximation to the vibrational degrees of freedom is not required in this new formulation, the vibrational-level dependencies of the transition frequencies and the vibrational relaxation rates are fully incorporated in the obtained formula for the 2D vibrational response function. Furthermore, the non-Markovian nature of the vibrational dephasing process in condensed phases is also fully taken into account by carrying out the resummation of the associated diagrams with the linked diagram theory. It is found that there exists an additional contribution to the 2D vibrational response function, which was completely ignored due to the approximations based on the Wick's theorem.