Ghost imaging via sparsity constraint (GISC) can recover objects from the intensity fluctuation of light fields even when the sampling rate is far below the Nyquist sampling rate. In this paper, we develop an efficient algorithm called the preconditioned multiple orthogonal least squares (PmOLS) for solving the GISC reconstruction problem. Our analysis shows that the PmOLS algorithm perfectly recovers any $n$-dimensional $K$-sparse signal from $m$ random linear samples of the signal with probability exceeding $1-3n^2 e^{-cm/K^2}$. Simulations and experiments demonstrate that the proposed algorithm has very competitive imaging quality compared to the state-ofthe-art methods.