Low-lying singlet excited states of pyrene derivatives originated from the ¹L_a and ¹L_b states of pyrene have decisive influences on their absorption and fluorescence emission behaviors. Calculation of these excited states with quantitative accuracy is required for the theoretical design of pyrene derivatives tailored to target applications; this has been a long-standing challenge for ab initio quantum chemical calculations. In this study, we explore an adequate computational scheme through calculations of pyrene and its phenyl-substituted derivatives using multi-reference perturbation theory (MRPT) methods. All valence π orbitals on the pyrene moiety were assigned to the active orbitals. Computational load was reduced by restricting the electron excitations within the active orbitals in the preparation of reference configuration space. A generalized multi-configuration quasi-degenerate perturbation theory (GMCQDPT) was adopted to treat the reference space other than the complete active space. The calculated ¹L_a and ¹L_b excitation energies of pyrene are in good agreement with the experimental values. Calculations of 1,3,6,8-tetraphenyl pyrene suggest that the energetic ordering of ¹L_a and ¹L_b is inverted through tetraphenyl substitution and its lowest singlet excited state is the ¹L_a parentage of pyrene, which is consistent with the experimentally deduced scheme. These results are not readily obtained by MRPT calculations with a limited number of active orbitals and single-reference theory calculations. Diphenyl pyrenes (DPPy) were also calculated at the same level of theory to investigate the dependence on the substitution positions of phenyl groups.