In this paper (paper I), we outline two general theoretical techniques. The first one is related to the sparsity of the HTC Hamiltonian, which allows one to reduce the cost for action of polariton Hamiltonian onto a state vector to the linear order in the number of states, instead of the quadratic order. The second one is applying the wellknown Chebyshev series expansion approach for quantum dynamics propagation and applying them 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 first technique we developed in this work. These two theoretical approaches are general and can be applied to any trajectory-based non-adiabatic quantum dynamics methods or gaussian wavepacket-based approaches. We focus on simulating the linear absorption spectra of the HTC model system, with both inhomogeneous site-energy disorder as well as dipolar orientational disorders. With the significant reduction of the computational costs, we can directly simulate collective polariton dynamics with up to N = 10 5 . Our numerical results agree well with the previous analytic 37 and numerical work.