The trick is that a translation of a point (a1, ..., an) is essentially a shear of a point (a1, ..., an, 1) - a linear operation in R^(n+1). But because the points that are of interest lie on a (hyper)plane that is not a vector space - it does not intersect the origin - a shear looks like a nonlinear operation after projection to R^n.

    ^
    |---a
    |   |
  --+--------->
    |

    ^
    | ,---a
    |/   /
  --+--------->
    |

And you can see how points on the a_(n+1) = 0 plane ("infinitely far away" in projective terms) are invariant under this operation - so they behave like pure direction vectors and can be used as normals and transformed with the same matrix as points!

This explanation really helped, thank you.

The translation part is not that complex, you're basically just introducing an extra coordinate 'w' which is always 1, so you can write x -> x + w, which is linear, rather than x -> x + 1, which is not.

Things start to get weird when you note that this works even when w isn't 1. If you then decide that e.g. (x,y,w) and (x/w, y/w, 1) should be the same point, you can do some rather interesting things. Although, if you identify (x:y:w) with all the points on the line (lx:ly:lw) (for any l) then suddenly (x/w, y/w, 1) is just the point where this intersects the plane w=1, which isn't too hard to visualize. And if you think about it this also kind of explains why this coordinate system can be used to project a 3D scene onto a 2D surface with correct perspective.