Over-the-air computation (AirComp) has been recognized as a promising technique in Internet-of-Things (IoT) networks for fast data aggregation from a large number of wireless devices. However, as the number of devices becomes large, the computational accuracy of AirComp would seriously degrade due to the vanishing signal-to-noise ratio (SNR). To address this issue, we exploit the massive multiple-input multiple-output (MIMO) with hybrid beamforming, in order to enhance the computational accuracy of AirComp in a cost-effective manner. In particular, we consider the scenario with a large number of multi-antenna devices simultaneously sending data to an access point (AP) equipped with massive antennas for functional computation over the air. Under this setup, we jointly optimize the transmit digital beamforming at the wireless devices and the receive hybrid beamforming at the AP, with the objective of minimizing the computational mean-squared error (MSE) subject to the individual transmit power constraints at the wireless devices. To solve the non-convex hybrid beamforming design optimization problem, we propose an alternating-optimization-based approach. In particular, we propose two computationally efficient algorithms to handle the challenging receive analog beamforming problem, by exploiting the techniques of successive convex approximation (SCA) and block coordinate descent (BCD), respectively. It is shown that for the special case with a fully-digital receiver at the AP, the achieved MSE of the massive MIMO AirComp system is inversely proportional to the number of receive antennas. Furthermore, numerical results show that the proposed hybrid beamforming design substantially enhances the computation MSE performance as compared to other benchmark schemes, while the SCA-based algorithm performs closely to the performance upper bound achieved by the fully-digital beamforming.