We consider the task of estimating the expectation value of an (n)-qubit tensor product observable (O_1\otimes O_2\otimes \cdots \otimes O_n) in the output state of a shallow quantum circuit. This task is a cornerstone of variational quantum algorithms for optimization, machine learning, and the simulation of quantum many-body systems. Here we study its computational complexity for constant-depth quantum circuits and three types of single-qubit observables (O_j) which are (a) close to the identity, (b) positive semidefinite, (c) arbitrary. It is shown that the mean value problem admits a classical approximation algorithm with runtime scaling as (\mathrm{poly}(n)) and (2^{\tilde{O}(\sqrt{n})}) in cases (a,b) respectively. In case (c) we give a linear-time algorithm for geometrically local circuits on a two-dimensional grid. The mean value is approximated with a small relative error in case (a), while in cases (b,c) we satisfy a less demanding additive error bound. The algorithms are based on (respectively) Barvinok’s polynomial interpolation method, a polynomial approximation for the OR function arising from quantum query complexity, and a Monte Carlo method combined with Matrix Product State techniques. We also prove a technical lemma characterizing a zero-free region for certain polynomials associated with a quantum circuit, which may be of independent interest.