Density functional theory plus Hubbard U (DFT+U) methodology was used to calculate the structures and energetic landscapes of CeSiO4, including its stetindite and scheelite phases from ambient pressure to ∼24 GPa. To ensure accurate simulations of the high-pressure structures, assessments of strain-stress methods and stress-strain methods were conducted in prior, with the former found to have a better agreement with the experimental result. From DFT calculations the equation of states (EOS) of both stetindite and scheelite were further obtained, with the fitted bulk moduli being 182(2) GPa and 190.0(12) GPa, respectively. These results were found to be consistent with the experimental values of 177(5) GPa and 222(40) GPa. Furthermore, the calculated energetics suggest that the stetindite structure is more thermodynamically stable than the scheelite structure at a pressure lower than 8.35 GPa. However, the stetindite → scheelite phase transition was observed experimentally at a much higher pressure of ∼15 GPa. A further phonon spectra investigation by the density functional perturbation theory (DFPT) indicated the Eg1 mode is being softened with pressure and becomes imaginary after 12 GPa, which is a sign of the lattice instability. Consequently, it was concluded that the stetindite → scheelite transition is predominantly initiated by the lattice instability under high-pressure.