By using stochastic geometry theorem, we investigate the spectral efficiency (SE) of largescale heterogeneous cellular networks (HetNets) with massive MIMO relays, where the maximum ratio combining/ maximal ratio transmission, zeroforcing processing scheme as well as amplifyandforward are employed at base stations and the macro cells are overlaid with dense small cells. To obtain the aggregate SE in target cell, we first formulate the signaltointerference noise ratios (SINR) for a user pair as well as the expression of power amplification factor. We then achieve the closedform expressions of aggregate SE. We find that when the radius of small cells is less than the guard radius we can acquire optimal (or suboptimal) system performance. For the guard radius, there is an upper bound, only under which the increase of the guard radius can contribute to the improvement of SE. Meanwhile, the radius of small cells has an optimal lower bound and only when the radius of macro cells is less than the optimal lower bound, the aggregate SE increases with the macro cell radius.