In this paper, we extended the L-PLDM formalism to simulate the non-linear spectra for N -molecules collectively coupled to the cavity, described by the HTC Hamiltonian. A direct density-matrix-based simulation will require a computational cost (in both time and memory requirements) proportional to O(N 6 ). The sparse nature of the HTC system Hamiltonian and dipole matrix let us compactly express ĤQ |Ψ s ⟩ and μ|Ψ s ⟩ as simple Hadamard products between vectors of O(N 2 ) instead of matrix-matrix or matrix-vector operations. We further utilize the forward-backward nature of the L-PLDM method to convert all matrix-matrix operations to consecutive vector-vector Hadamard products. Here, we incorporated these strategies into the previously developed L-PLDM approach for simulating 2D Electronic Spectroscopy (2DES) spectra of exciton-polariton under the collective coupling regime. In particular, we apply the efficient quantum dynamics propagation scheme developed in Paper I 1 to both the forward and the backward propa- gations in the PLDM, and develop an efficient important sampling scheme and GPU vectorization scheme that allow us to systematically reduce the computational costs from O(K 2 )O(T 3 ) to O(K)O(T 0 ) for the 2DES spectra simulation, where K is the number of states and T is the number of time steps of propagation. We also simulate the 2DES of N = 25 molecules coupled to the cavity within the double excitation subspace and demonstrate the effect of linewidth reduction due to the polaron decoupling effect, which has been observed experimentally. We further analyzed the signal from both rephasing and non-rephasing contributions, as well as the ground state bleaching (GSB), excited state emission (ESA), and stimulated emission (SE). For N > 1, we also found the derivative lineshape due to ESA, which appears as a negative feature on the lower and upper polariton diagonal peaks due to the difference in frequencies of |G⟩ → |P 1 ⟩ and |P 1 ⟩ → |P 2 ⟩ transitions. Together with paper I, the theoretical development we presented here allows an efficient and accurate simulation of quantum dynamics and multi-dimensional spec-troscopy of exciton polariton under the collective coupling regimes, which we envision will benefit the investigation of molecular polariton photophysics.