Another issue is that the system is overdetermined (four equation in three variables) and has no generic solution. When I remove the radicals by squaring I get reasonable speed and an empty solution set.
Solve[{rps^2 == (xs - xp)^2 + (ys - yp)^2 + (zs -
zp)^2, (rps Sin[theta] Cos[phi] + xp)^2 ==
4 f zs + 4 f^2 - ys^2, (rps Sin[theta] Sin[phi] + yp)^2 ==
4 f zs + 4 f^2 - xs^2,
rps Cos[theta] + zp == (xs^2 + ys^2)/(4 f) - f}, {xs, ys, zs}]
(* Out[385]= {} *)
