This paper considers pilot-based channel estimation in large-scale multiple-input multiple-output (MIMO) communication systems, also known as "massive MIMO". Unlike previous works on this topic, which mainly considered the impact of inter-cell disturbance due to pilot reuse (so-called pilot contamination), we are concerned with the computational complexity. The conventional minimum mean square error (MMSE) and minimum variance unbiased (MVU) channel estimators rely on inverting covariance matrices, which has cubic complexity in the multiplication of number of antennas at each side. Since this is extremely expensive when there are hundreds of antennas, we propose to approximate the inversion by an L-order matrix polynomial. A set of low-complexity Bayesian channel estimators, coined Polynomial ExpAnsion CHannel (PEACH) estimators, are introduced. The coefficients of the polynomials are optimized to yield small mean square error (MSE). We show numerically that near-optimal performance is achieved with low polynomial orders. In practice, the order L can be selected to balance between complexity and MSE. Interestingly, pilot contamination is beneficial to the PEACH estimators in the sense that smaller L can be used to achieve near-optimal MSEs.