Homogeneity of space implies linearity of Lorentz transformations
I claim that if the transformation between frames is homogeneous and differentiable, then it is affine (homogeneity is not strictly speaking sufficient for linearity since the full transformation between frames is actually a Poincare transformation which is affine, not linear)
For a mathematically precise proof, we need a mathematical definition of homogeneity. To arrive at such a definition, we note that the basic idea is that we can pick our origin wherever we choose, and it won't "affect the measurement results of different observers." In particular, this applies to measurements of the differences between the coordinates of two events. Let's put this in mathematical terms.
Let $L:\mathbb R^4\to\mathbb R^4$ be a transformation. We say that $L$ is homogeneous provided \begin{align} L(x+\epsilon) - L(y+\epsilon) = L(x) - L(y) \end{align} for all $\epsilon\in\mathbb R^4$ and for all $x,y\in\mathbb R^4$.
We can now precisely state and prove the desired result. Note that I also assume that the transformation is differentiable. I haven't thought very hard about if or how one can weaken and/or motivate this assumption.
Proposition. If $L$ is homogeneous and differentiable, then $L$ is affine.
Proof. The definition of homogeneity implies that , \begin{align} L(x+\epsilon)-L(x) = L(y+\epsilon) - L(y) \tag{1} \end{align} for all $\epsilon, x, y$. Now we note that the derivative $L'(x)$ of $L$ at a point $x$ is a linear operator on $\mathbb R^4$ that satisfies \begin{align} L(x+\epsilon) - L(x) = L'(x)\cdot \epsilon +o(|\epsilon|) \end{align} and plugging this into $(1)$ gives \begin{align} (L'(x)-L'(y))\cdot\epsilon = o(|\epsilon|) \end{align} for all $\epsilon,x,y$, where $|\cdot|$ is the euclidean norm. Now simply choose $\epsilon = |\epsilon|e_j$ with $|\epsilon|\neq 0$ where $e_0, \dots e_3$ are the standard, ordered basis elements on $\mathbb R^4$, multiply both sides on the left by $(e_i)^t$ where $^t$ means transpose, divide both sides by $|\epsilon|$, and take the limit $|\epsilon|\to 0$ to show that all matrix elements of $L'(x)-L'(y)$ are zero. If follows immediately that \begin{align} L'(x) = L'(y) \end{align} In other words, the derivative of $L$ is constant. It follows pretty much immediately that $L$ is affine, namely that there exists a linear operator $\Lambda$ on $\mathbb R^4$, and a vector $a\in\mathbb R^4$ such that \begin{align} L(x) = \Lambda x + a \end{align} for all $x\in\mathbb R^4$. $\blacksquare$
This answer is essentially the same as JoshPhysics's Answer but with the following points:
- We use a more "off-the-shelf" mathematical result to get rid of the differentiability assumption JoshPhysics used in his answer and instead we simply need to assume that the Lorentz transformation is only continuous;
- We show that the continuity assumption is a necessary and the minimum necessary assumption further to the OP's homogeneity assumption i.e. the OP's assertion that linearity comes from the homogeneity of spacetime alone is wrong.
JoshPhysics's equation (1) implies:
$$L(X+Y)-L(Y) = L(X)-L(0);\;\forall X,\,Y\in \mathbb{R}^{1+3}\tag{1}$$
Now we define $h:\mathbb{R}^{1+3}\to\mathbb{R}^{1+3}$ by $h(Z)=L(Z)-L(0)$; then it follows from (1) alone that:
$$h(X+Y)=h(X)+h(Y);\quad\forall\,X,\,Y\in\mathbb{R}^{1+3}\tag{2}$$
But this is the famous Cauchy functional equation generalized to $3+1$ dimensions. For one, real dimension, the only continuous solution is $h(X)\propto X$; there are other solutions, but they are everywhere discontinuous, as shown in:
E. Hewitt & K. R. Stromberg, "Real and Abstract Analysis" (Graduate Texts in Mathematics), Springer-Verlag, Berlin, 1965. Chapter 1, section 5
It is easy to broaden the Hewitt-Stromberg argument to any number of dimensions, so that, given an assumption of continuity of $L:\mathbb{R}^{1+3}\to\mathbb{R}^{1+3}$, we must have:
$$L(X) = \Lambda\,X + \Delta\tag{3}$$
where $\Lambda$ is a linear operator - a $4\times4$ matrix and $\Delta$ a spacetime offset.
Note that we must invoke the continuity assumption; otherwise, following the reasoning in Hewitt and Stromberg, our $h$ function could be one of the everywhere discontinuous Cauchy equation solutions, and we could then, by reversing the step from my equation (1) to (2), construct everywhere discontinuous, nonlinear functions that fulfill JoshPhysics's homogeneity postulate. So unless we require exactly continuity, we shan't "choose" the right solution of the Cauchy equation. Thus continuity of transformation as well as homogeneity are the minimum assumptions needed to imply linearity.