In this paper, we consider the problem of channel estimation for large intelligent metasurface (LIM) assisted massive multiple-input multiple-output (MIMO) systems. The main challenge of this problem is that the LIM with a large number of low-cost metamaterial antennas can only passively reflect the incident signal by a certain phase shift, and does not have any signal processing capability. We introduce a general framework for the estimation of the transmitter-LIM and LIM-receiver cascaded channel, and propose a three-stage algorithm that includes a sparse matrix factorization stage, a probabilistic ambiguity elimination stage, and a matrix completion stage. Simulation results illustrate that the proposed algorithm can achieve accurate channel estimation for LIM-assisted massive MIMO systems.