In recent years, phase retrieval has received much attention in statistics, applied mathematics and optical engineering. In this paper, we propose an efficient algorithm, termed Subspace Phase Retrieval (SPR), which can accurately recover an $n$-dimensional $k$-sparse complex-valued signal $\x$ given its $Ω(k^2\log n)$ magnitude-only Gaussian samples if the minimum nonzero entry of $\x$ satisfies $|x_{\min}| = Ω(\|\x\|/\sqrt{k})$. Furthermore, if the energy sum of the most significant $\sqrt{k}$ elements in $\x$ is comparable to $\|\x\|^2$, the SPR algorithm can exactly recover $\x$ with $Ω(k \log n)$ magnitude-only samples, which attains the information-theoretic sampling complexity for sparse phase retrieval. Numerical Experiments demonstrate that the proposed algorithm achieves the state-of-the-art reconstruction performance compared to existing ones.