In this paper, we presented two general theoretical techniques for efficient simulation of polariton quantum dynamics under the collective light-matter coupling regime. The first one is related to the sparsity of the HTC Hamiltonian, which allows one to reduce the cost of acting polariton Hamiltonian onto a state vector to the linear order of the number of states, instead of the quadratic order. In principle, a direct density-matrix-based simulation requires a computational cost (in both time propagation and memory requirements) proportional to O(N 6 ), where N is the number of molecules and O(N 2 ) is the number of states in the double excitation subspace. The sparse nature of the HTC system Hamiltonian and dipole matrix allow for a compact expression of ĤQ |Ψ s ⟩ and μ|Ψ s ⟩ as simple Hadamard products between vectors of O(N 2 ) instead of matrix-matrix or matrix-vector operations.