For the last several decades, there has been a rapid development of nanoscience and nanotechnology. In particular, nanoparticles (NPs) are applied in various catalyst problems, where their enormously large surface-to-volume ratio not only is advantageous for high catalytic performance but also allows the quantum size effect to play a key role in modifying their chemical properties. However, when understanding the size effect of NPs on the catalytic properties by employing density functional theory (DFT) calculations, there has been an obvious experiment-theory gap in simulating nanocatalysts with realistic sizes. In this study, we developed a new simulation method based on the cluster expansion model, namely, CE-np, which enables efficient and accurate calculations of the intermediate binding energies for various sizes of NPs. We then applied CE-np to investigate the electrochemical CO2 reduction reaction (CO2RR) of gold NPs (AuNPs). CE-np reproduces not only DFT-level accuracies in predicting the intermediate binding energies on the NPs and slab surfaces but also the experimental behavior of catalytic activity and selectivity of AuNP catalysts. Because of the high computational efficiency of CE-np (without sacrificing the accuracy level of DFT), we performed the most exhaustive search on all possible on-top binding sites of AuNPs to unveil the complicated relations between the catalytic performance (activity and selectivity) and the NP properties (shape and size). This also highlights for the first time the catalytic importance of the near-edge sites, that is, active sites on the facets that are very close to the edges. We anticipate that our methodological development and several new findings on the CO2RR activity of AuNPs will provide advances in developing CO2 electrochemical reduction technologies based on NPs.