We develop two iterative algorithms for solving the low rank phase retrieval (LRPR) problem. LRPR refers to recovering a low-rank matrix $\X$ from magnitude-only (phaseless) measurements of random linear projections of its columns. Both methods consist of a spectral initialization step followed by an iterative algorithm to maximize the observed data likelihood. We obtain sample complexity bounds for our proposed initialization approach to provide a good approximation of the true $\X$. When the rank is low enough, these bounds are significantly lower than what existing single vector phase retrieval algorithms need. Via extensive experiments, we show that the same is also true for the proposed complete algorithms.