Large-scale variational quantum algorithms are widely regarded as a potential pathway to achieve practical quantum advantages. However, the presence of quantum noise might suppress and undermine these advantages, which blurs the boundaries of classical simulability. To gain further clarity on this matter, we present a novel polynomial-scale method based on path integrals on the Pauli basis. This method efficiently approximates quantum mean values in variational quantum algorithms with bounded truncation error in the presence of independent single-qubit depolarizing noise. Theoretically, we have rigorously proved: 1) For a fixed noise rate $\lambda$, our method's time and space complexity exhibit a polynomial relationship with the number of qubits $n$, the circuit depth $L$, the inverse truncation error $\frac{1}{\varepsilon}$, and the root square inverse success probability $\frac{1}{\sqrt{\delta}}$. 2) For variable $\lambda$, the computational complexity becomes $\mathrm{Poly}\left(n,L\right)$ when $\lambda$ exceeds $\frac{1}{\log{L}}$ and it becomes exponential with $L$ when $\lambda$ falls below $\frac{1}{L}$. Numerically, we have conducted classical simulations of IBM's zero-noise extrapolated experimental results on 127-qubit Eagle processor [Nature \textbf{618}, 500 (2023)]. Our approach attains higher accuracy and faster runtime compared to the quantum device. Moreover, this method enables us to deduce noisy outcomes from noiseless results, allowing us to accurately reproduce IBM's unmitigated results that directly correspond to experimental observations.