Theorem Schmeorem. A Galilean invariant Lagrangian for any number of classical particles interacting with a potential:
$$ S = \int \sum_k {m_k(\dot{x}_k-u)^2\over 2} + \lambda \dot{u} - U(x_k)\;\;\; dt $$
For any Galilean invariant Lagrangian $L(\dot{x}_k, x_k)$, the Lagrangian
$$ L'(\dot{x}_k,x_k, \lambda, u) = L(\dot{x}_k-u,x_k) + \lambda \dot{u} $$
is explicitly Galilean invariant, and has the same dynamics (assuming the original Lagrangian was Galilean invariant).
The Galilean properties of the x's are as usual. The dynamical variables extend to include $\lambda,u$ which act as Lagrange multipliers. The transformation law for u and $\lambda$ are:
$ x \rightarrow x-vt $
$ u \rightarrow u-v $
$ \lambda \rightarrow \lambda $
And it is trivial to verify that the new Lagrangian is completely invariant. The equation of motion for $\lambda$ just makes $u$ constant, equal to $u_0$, while the equation of motion for $u$ integrates to
$$ \lambda = - \sum_k m_k x_k - M u_0 t $$
up to an additive constant which I have set to zero. This is almost all the equations of motion, but there is one more equation which comes from extremizing the action with respect to $u_0$, which sets
$$ u_0 = \sum_k m_k \dot{x}_k $$
Where the time is unimportant, because this is the center of mass velocity, which is conserved. The Noether prescription in the explicitly Galilean invariant action is trivial--- the conserved quantity associated with Galilean boosts is just $\lambda$, and this is indeed the center of mass position.
Why this works
If you integrate the kinetic energy for the usual free particle action by parts, you get:
$$S = \int \sum_k m\ddot{x}_k x_k + U(x_k) dt$$
This action is Galilean invariant on mass shell, meaning that the non-galilean invariant part is zero when you enforce the equations of motion. This means that adding some additional nondynamical fields should produce a Galilean invariant action off shell, and this is the $\lambda, u$.
Relation to Lorentz transformations
When you perform a Lorentz transformation, the arclength particle action is invariant. But if you fix the origin of the Lorentz transformation on the initial time, the final time is transformed, so the path is not going to the same final time anymore after the transformation. When you take the non-relativistic limit, the final time becomes degenerate with the initial time, but the action cost from shifting the final time does not approach zero.
This means that you need an extra variable to keep track of the infinitesimal bit of final time, and that this extra variable will need a nontrivial transformation law under Galilean transformations.
To find out what this new variable should be, it is always best to consider the analogous thing for rotational invariance. Consider a string at tension with small deviations from horizontal, and let the deviation of the string from horizontal be h(t). The rotationally invariant potential energy is the arclength of the string
$$ U = \int \sqrt{1+h'^2} dx $$
and this is the potential energy which gives the rotationally invariant analog of the wave equation. Once you go to small deviations, the expansion for U gives the usual wave-equation potential energy
$$ U(h) = \int {1\over 2} h'^2 dx $$
and this is no longer rotationally invariant. But it is skew-invariant, meaning adding a constant slope line to h does not change the energy. Except that it does, by a perfect derivative:
$$ U(h + ax) = \int {1\over 2}h'^2 + a h' + {a^2\over 2} dx$$
This is clearly the same exact situation as for the Lorentz invariance turning into Galilean invariance, except using rotational invariance, where everybody's intuition is firm. The additional $a^2\over 2$ energy is due to the quadratic extra length of a rotated string, while the linear perfect derivative $ah'$ integrates to $a (h_f - h_i)$, and this is the amount of reduction/increase in length when you rotate a tilted string.
So to get a fully tilt-invariant potential energy, you need to add a variable $u$ which is dynamically constrained to equal the total tilt of the string. This variable will distinguish between different rotated versions of the string: rotating the string by itself without rotating the average tilt variable will change the energy--- this is because tilting the horizontal string between 0 and A is not quite the same as the pre-tilted string between 0 and A, the pre-tilted string has a different length. Rotating the total tilt by itself will change the energy, but rotating both does nothing, and this is the encoding of rotatonal invariance.
So you need an average tilt variable to turn explicit rotational invariance to explicit tilt invariance. The total potential energy is then given by the deviations from the average tilt:
$$ U = \int {1\over 2} (h'-u)^2 dx $$
and u transforms as $u-a$ under a tilt by a. This makes the potential energy invariant.
The kinetic energy is given by the time dependence of h, and there must be a Lagrange multiplier to enforce that the total tilt is equal to the average tilt
$$ S = \int {1\over 2} \dot h^2 - {1\over 2} (h'-u)^2 + \beta (u - h') dt dx $$
Where $\beta$ is a global in x Lagrange multiplier for u, forcing it to equal h'. But it does no harm to allow u to vary in x, so long as the Lagrange multiplier enforces that it is constant. The way to do this is to change the Lagrange multiplier term to
$$ - \int \lambda' (u(x) - h'(x)) dx = \int \lambda (u'(x) - h''(x)) $$
But then the equation of motion kills the second term, so you only need a Lagrange multiplier to be:
$$ \int \lambda u'(x)$$
And the equations of motion automatically constrain u to be the average slope. These manipulations have exact analogs in Lorentz transformations, and explain the relation of the explicitly Galilean invariant action to the Lorentz action. The analog of average slope is the center of mass velocity.
Best Answer
Yes, the invariance of the action follows from special relativity – and special relativity is right (not only) because it is experimentally verified. All the equations of motion may be derived from the condition $\delta S = 0$, the action is stationary (which usually means it has the minimum value on the allowed trajectory/history among all trajectories/histories with the same initial and final conditions). If $S$ depended on the inertial system, so would the terms in the equations $\delta S =0$, and these laws of motion couldn't be Lorentz-covariant (note how this Lorentz is spelled; Lorenz also existed but it was a different physicist). Quite generally, you shouldn't think about "derivation of the action". When we work with the action at all, we are doing so because we view the action as the most fundamental expression – and we derive everything else out of it. In that context, we pretty much define a Lorentz-invariant theory as a theory determined by a Lorentz-invariant action.
Your integral is Lorentz-invariant but it is not translationally invariant under $x^\mu \to x^\mu + a^\mu$. So it's not Poincaré-invariant (the Poincaré symmetry unifies the Lorentz transformations and spacetime translations) and due to this violation, we also say that it disagrees with the laws of special relativity. You could also create other expressions, e.g. replace $x_\mu x^\mu$ in the integral by some extrinsic curvature invariant of the world line etc. Those terms could be made Poincaré-invariant. So the right claim is that the proper length of the world line is the only Poincaré-invariant functional that doesn't depend on any higher derivatives of the coordinates $x^\mu(\tau)$.