The reduced density matrix (RDM) is crucial in quantum many-body systems for understanding physical properties, including all local physical quantity information. This study aims to minimize various error constraints that causes challenges in higher-order RDMs estimation in quantum computing. We identify the optimal balance between statistical and systematic errors in higher-order RDM estimation in particular when cumulant expansion is used to suppress the sample complexity. Furthermore, we show via numerical demonstration of quantum subspace methods for one and two dimensional Fermi Hubbard model that, biased yet efficient estimations better suppress hardware noise in excited state calculations. Our work paves a path towards cost-efficient practical quantum computing that in reality is constrained by multiple aspects of errors.