We demonstrate strategies to efficiently manipulate matrix-vector multiplications to reduce the computation cost of semiclassical simulations. When considering the single excitation subspace, the size of the Hilbert space K ∝ O(N ), and solving the TDSE in principle requires operation of O(K 2 ) ∝ O(N 2 ). When considering the Hilbert space up to the second excitation subspace (for 2DES simulations in paper II) K ∝ O(N 2 ), and the quantum dynamics simulation will in principle require a cost of O(K 2 ) ∝ O(N 4 ). To reduce the scaling of the cost, we first demonstrate how the action of the ĤQ on |Ψ(t)⟩ can be written as the sum of a simple Hadamard product of vectors that reduces matrix-vector multiplication from O(N 2 ) to O(N ) for single excitation subspace and O(N 4 ) to O(N 2 ) for up to double excitation subspace. We demonstrate a similar scaling reduction for the action of the dipole operator on the wavefunction, which is relevant for linear and non-linear spectra calculations. We further apply a Chebyshev series expansion approach to simulate the polariton dynamics in the HTC system. The Chebyshev expansion allows us to use a much larger time step for propagation and only requires a few recursive operations of the Polariton Hamiltonian acting on state vectors, which can further take advantage of the scaling reduction due to the sparsity of the HTC Hamiltonian.