This paper investigates the transmitter and receiver beamforming (TB and RB) for over-the-air computation (AirComp) in massive multiple-inputs and multiple-outputs (MIMO) systems. First, we propose a two-phase hybrid beamforming algorithm to design TB and hybrid RB. In the first phase, we adopt a projected gradient descent with momentum (PGDM) algorithm to search for the optimal fully-digital TB matrices. Compared with the benchmarks on the mean square error (MSE) performances, PGDM can achieve up to 5 dB gain in signal-to-noise ratio (SNR) with less algorithm execution time when fully-digital RB is assumed. In the second phase, we plug the TB matrices obtained in PGDM as well as the optimal baseband RB (BBRB) matrix into the MSE objective, and adopt gradient descent to search for the optimal radio-frequency RB (RFRB) matrix. Compared with the state-of-the-arts, the proposed two-phase algorithm reduces up to 30% of the algorithm execution time and 18% of the MSE. Second, we propose a statistical TB algorithm to reduce the communication overheads, in which TB completely depends on statistical channel state information (CSI) and thus does not rely on the feedback from the base station (BS). We prove that orthonormal matrices are asymptotically optimal for statistical TB when uncorrelated Rayleigh channels are assumed and the number of receiving antennas approaches to infinity. For correlated channels, experimental results show that the proposed statistical TB can achieve about 5 dB gain in SNR compared with orthonormal matrices in terms of the MSE performance. Third, a large scale system analysis is made in this paper. As the number of the receiving antennas approaches to infinity, asymptotically optimal choices for TB and RB are provided, and upper bounds of MSE are derived in terms of the number of clients and receiving antennas. For hybrid RB, an upper bound of the squared distance between the optimal hybrid RB and the optimal fully-digital RB is also derived.