"A not so simple approach". I may have messed a little bit too much with grouping the terms, do forgive my elementary math skills, it's a side effect of using tools like wolfram and mathematica too much.
Without loss of generality we can assume that the near plane is given as $z=1$ (that is having normal $(0,0,1)$ and a point on it $(0,0,1)$) and the camera center (the reflected pinhole camera aperture point) is at $(0,0,0)$. Let $p,q \in \mathbb{R}^3$ be two points on a line (non-degenerate, meaning $p\ne q$, otherwise it's just a point and the projection is also a point - unless it's at $(0,0,0)$ then it's undefined). If the line defined by $p,q$ passes through $(0,0,0)$ the projection is a point and the requirement is trivially satisfied.
Now let this not be the case. A perspective projection onto the plane through $(0,0,0)$ is given by $p' = \frac{p}{p_z}$, and $q' = \frac{q}{q_z}$. The line defined by $p,q$ is given as $r(\lambda) = p + \lambda(q-p), \lambda \in \mathbb{R}$. Its projection is given by $r'(\lambda) = \frac{r(\lambda)}{r_z(\lambda)}$. Since $z=1$ for all projected points, then we need only study the $x,y$ coordinates. We want to show that any projected point $r'(\lambda)$ can be written as an affine combination of $p',q'$ (this defines a line). Taking that into account:
$$p'+\mu(q'-p') = \frac{r(\lambda)}{r_z(\lambda)}$$
$$(1-\mu)\frac{p}{p_z} + \mu\frac{q}{q_z} = \frac{p + \lambda(q-p)}{p_z + \lambda(q_z-p_z)}$$
$$(1-\mu)q_z(p_z+\lambda(q_z-p_z))p + \mu p_z(p_z+\lambda(q_z-p_z))q = q_zp_z(p + \lambda(q-p))$$
$$q_zp_z(p - \lambda p) + q^2_z\lambda p - \mu q_z(p_z+\lambda(q_z-p_z))p + \mu p_z(p_z+\lambda(q_z-p_z))q - q_zp_z(p - \lambda p) - q_zp_z\lambda q = 0$$
$$(1-\mu)\lambda q^2_z p -\mu(1-\lambda)q_zp_z p + \mu(1-\lambda)p^2_zq - (1-\mu)\lambda p_zq_z q = 0$$
$$(1-\mu)\lambda q_z(q_zp - p_zq) - \mu(1-\lambda)p_z(q_zp-p_zq) = 0$$
$$((1-\mu)\lambda q_z - \mu(1-\lambda)p_z)(q_zp - p_zq) = 0$$
For this to be satisfied one of the terms must be $0$. The first case concerns the second term: $p' = q'$, which we already considered as the trivial case. Then we must show that we can always pick a $\mu$ such that $(1-\mu)\lambda q_z - \mu(1-\lambda)p_z = 0$. Solve for $\mu$ and you get: $\mu = \frac{\lambda q_z}{(1-\lambda)p_z + \lambda q_z}$. Thus I have not only found that it is a straight line, but also the relationship between the parameterizations.
In general a projection onto a different 2-manifold (not a plane) will not result in the lines remaining straight (think of a sphere or a cylinder).
