I have removed my earlier answer which relied on Taylor expansions, and will try and give the argument without them (they are of course equivalent) including how these arguments extend to general relativity and comparisons to the Euclidean case. If this makes sense I can restore the alternative Taylor expansion argument, but the references and video at the end also contain it.
First, let's examine some of the Euclidean versions of the SR arguments to make sure we understand what's going on. In Euclidean space we know that
$$ds = \sqrt{dx^2 + dy^2 + dz^2}$$
is just the length of a(n infinitesimal) line segment in 3D space, a real number that we can treat as the argument of a function or as a function of some other variable, or treat it as a parameter/the-variable in a vector function etc... - it's just the Pythagorean theorem for a small line segment. We can even introduce a metric tensor $g_{ij}$ by writing
$$ds = \sqrt{dx^2 + dy^2 + dz^2} = \sqrt{d \mathbf{r} \cdot d \mathbf{r}} = \sqrt{\frac{\partial \mathbf{r}}{\partial u^i} \cdot \frac{\partial \mathbf{r}}{\partial u^j} du^i du^j} = \sqrt{g_{ij} du^i du^j}$$
but the interpretation of $ds$ as a simple line segment has gone nowhere. It barely matters that this involves a square root of a quadratic form in a vector space, it's just the Pythagorean theorem nothing else. Since vector spaces assume every vector is situated at the origin, but we are clearly thinking of $d \mathbf{r} = (dx,dy,dz)$ as the tangent vector to a curve at an arbitrary point along the curve, it's going to involve a lot of ultimately trivial technicalities all just so we can (redundantly) attach a vector space to each point of another 'vector space' (=space-time, by interpreting this vector space as a manifold, because we don't care about adding the vectors representing the points of space-time, it makes it look less silly than placing a vector space at each point of a vector space) and to then define a metric in each separate vector space just to be able to recreate the Pythagorean theorem. We could just talk in terms of affine spaces, but lets ignore this too, it all really just involves basically nothing more than the Pythagorean theorem and the thinking of setting up single-variable arc length integrals.
Here time $t$ is just a parameter allowing for a parametrization of $(x,y,z)$, it's a parameter like $s$ (or $ds$) is a parameter. In either case, $(x(t),y(t),z(t))$ or $(x(s),y(s),z(s))$, the result is a curve, and I can write
$$ds = \frac{ds}{dt} dt = \sqrt{\mathbf{v} \cdot \mathbf{v}} dt = \sqrt{g_{ij} \dot{u}^i \dot{u}^j} dt.$$
If I denote $\mathbf{v} \cdot \mathbf{v}$ as $a(v)$, and square the last equation, I have
$$ds^2 = a(v) dt^2.$$
We're going to see this notation below. The point is, if $s$ is a parameter whose differential is given in the first equation above, and $t$ is some other parameter (which allows us to interpret $\frac{d \mathbf{r}}{dt}$ as a velocity), then the relationship between the squares of their differentials in this case is mediated by a function of the magnitude of the velocity $v = \sqrt{\mathbf{v} \cdot \mathbf{v}}$ (we could obviously also denote the dependence on $(x,y,z)$ in $a(v)$).
I can obviously work with $ds^2$ instead of $ds$, where $ds^2$ is also just a parameter, it doesn't matter that it's also given as a quadratic form relative to some coordinate system, I am interpreting it as a single parameter, a real number. If I write
$$ds^2 = dx^2 + dy^2 + dz^2 = [(\frac{dx}{dt})^2 + (\frac{dy}{dt})^2 + (\frac{dz}{dt})^2]dt^2 = v^2 dt^2$$
it becomes obvious I am just treating it as a parameter, that's what the $ds^2$ notation signifies. This is extremely standard notation for (the square of) an arc length segment found even in basic calculus books and that is clear from the context of L&L's discussion.
I can also obviously work with single-variable functions of $ds$, i.e. $f = f(ds)$, or of $ds^2$, i.e. $f = f(ds^2)= f(v^2 dt^2)$. Any-time $ds^2$ is written below or in L&L, one should understand it as the square of an infinitesimal distance between two points, a scalar quantity in the terminology of tensors, but we're not assuming anything about transformation laws and invariance's yet, for good reason as we'll see. If we Taylor expand in $ds$ or $ds^2$ we are not Taylor expanding in 'quadratic forms' or something, we are Taylor expanding in a single scalar parameter just like we would if we Taylor expanded a function of time $f = f(t) = f(t - 0) = f(dt)$. The Padmanabhan references below in particular frame the argument in terms of a Taylor expansion, as did my earlier answer, but below I will try to ignore this completely.
I could now continue the Euclidean discussion and do here what I'm about to do for SR and GR but, because time $t$ is now just an external parameter, it's more geometric so I will sketch it at the end of the post.
So let's now consider the relativistic case, both for special relativity (SR) and general relativity (GR) to really make the argument crystal clear. First, SR:
The principle of relativity says that the laws of physics are the same in all inertial frames. Let us define the line element in an inertial frame $K$ as
$$ds_I^2 = c^2 dt^2 - dx^2 - dy^2 - dz^2.$$
I can also define the quantity $ds_I'^2$ in some other (not necessarily inertial) frame $K'$ as
$$ds_I'^2 = c^2 dt'^2 - dx'^2 - dy'^2 - dz'^2$$
So far, the only thing I can say about the parameter $s_I$ in $K$ is that $ds_I^2 = 0$ for two points that are connected by a light beam as viewed from $K$. Indeed this was the motivation to work with $ds_I^2$ in the first place. If $K'$ is also inertial then $ds_I'^2 = 0$ also holds for two points that are connected by the same light beam as viewed from $K$, however if $K'$ is non-inertial we will see this need not be the case. I remind you that e.g. $x_2 = (ct,0,0,0)$ and $x_1 = (0,0,0,0)$ in $K$ are such that $ds_I^2 = c^2 dt^2 \neq 0$ so not every two points are connected by a light beam in $K$.
I can express $ds_I^2$ in terms of an arbitrary parameter $u'$, whose form I leave completely general for the moment, as
\begin{align}
ds_I^2 &= c^2 dt^2 - dx^2 - dy^2 - dz^2 \\
&= [(c \frac{dt}{du'})^2 - (\frac{dx}{du'})^2 - (\frac{dy}{du'})^2 - (\frac{dz}{du'})^2]du'^2 \\
&= (\frac{ds_I}{du'})^2 du'^2.
\end{align}
Given that we want to think of $ds_I^2$ as the analogue of the 3D $ds^2$ squared line element in $K$, i.e. an invariant quantity representing the distance between two infinitesimally close points (which we can integrate to then generate a curve), we would like it to be able to retain this interpretation in a new frame $K'$, where after expressing $ds_I^2$ in the $K'$ frame in terms of the $(t',x',y',z')$ coordinates by using to-be-determined (in general, non-linear) functions $t = t(t',x',y',z')$ etc... and then treating the new $(t',x',y',z')$ as a function of the parameter $u'$ in this new $K'$ frame. The question is, what parameter $u' = \int du'$ in the (potentially non-inertial) frame $K'$ is the right parameter which allows us to say that $ds_I^2 = du'^2$. Everything so far still applies even to non-inertial $K'$ frames (with $ds_I'^2$ still defined as above). We thus examine what is the meaning of this $(\frac{ds_I}{du'})^2$ term.
First, we note that the RHS will have to vanish when the LHS vanishes, but when $ds_I^2 = 0$ we see that $(\frac{ds_I}{du'})^2 du'^2 = 0$ could be due to either $(\frac{ds_I}{du'})^2$ or $du'^2$. To simplify things, the first thing we do is try to choose the parameter $u'$ to be such that $du'$ itself vanishes whenever $ds_I$ vanishes.
Only now do I need to interpret $K'$ as an inertial frame to continue my argument. I am still allowing for the possibility of non-linear transformations in the inertial case, and later we will see what happens if $K'$ is non-inertial.
Now that we assume $K'$ is inertial, we can invoke the principle of invariance of the speed of light in inertial frames to say that since $ds_I'^2 = c^2 dt'^2 - dx'^2 - dy'^2 - dz'^2$ vanishes whenever $ds_I^2 = c^2 dt^2 - dx^2 - dy^2 - dz^2$ vanishes, if $du'^2$ is a function of $ds_I'^2$ then we can reduce the vanishing of the RHS (when the LHS vanishes) to the vanishing of $du'^2 = du'^2(ds_I'^2)$ by choosing $u'$ to be such that
$$du'^2 = \frac{du'^2}{ds_I'^2} ds_I'^2$$
and re-write the above as
\begin{align}
ds_I^2 &= c^2 dt^2 - dx^2 - dy^2 - dz^2 \\
&= (\frac{ds_I}{du'})^2 du'^2 = (\frac{ds_I}{ds_I'})^2 (\frac{ds_I'}{du'})^2 (\frac{du'}{ds_I'})^2 ds_I'^2 = (\frac{ds_I}{ds_I'})^2 ds_I'^2 \\
&= [(c \frac{dt}{ds_I'})^2 - (\frac{dx}{ds_I'})^2 - (\frac{dy}{ds_I'})^2 - (\frac{dz}{ds_I'})^2]ds_I'^2 \\
&= (\frac{ds_I}{ds_I'})^2 ds_I'^2.
\end{align}
We now examine $(\frac{ds_I}{ds_I'})^2$. We can immediately see from
$$(\frac{ds_I}{ds_I'})^2 = (c \frac{dt}{ds_I'})^2[1 - (\frac{1}{c}\frac{dx}{dt})^2 - (\frac{1}{c}\frac{dy}{dt})^2 - (\frac{1}{c}\frac{dz}{dt})^2] = (c \frac{dt}{ds_I'})^2[1 - (\frac{v}{c})^2]$$
that this factor now seems to depend on velocity in one frame. It could also abstractly depend on other things like the coordinates $(t,x,y,z)$ and the velocities (we'll see this in the non-inertial example below). It is only due to the fact that both frames are inertial frames that we can now also say that no, it can't depend on anything other than the magnitude of the relative velocity between the frames due to isotropy of space and homogeneity of space and time (i.e. the symmetries underlying the usual additive Noether symmetries). Compare to the non-inertial example below. So we can define this factor as
$$a(v) = (\frac{ds_I}{ds_I'})^2 = (c \frac{dt}{ds_I'})^2[1 - (\frac{1}{c}\frac{dx}{dt})^2 - (\frac{1}{c}\frac{dy}{dt})^2 - (\frac{1}{c}\frac{dz}{dt})^2] = (c \frac{dt}{ds_I'})^2[1 - (\frac{v}{c})^2]$$
and now say that
$$ds_I^2 = a(v) ds_I'^2.$$
It's still not certain that $ds_I'^2$ is the right parameter, only if $a(v) = 1$ is it the right parameter. That is what the rest of L&L's argument is: since $a(v)$ only depends on the magnitudes of the relative velocities, it is natural to examine the relative velocities between two different frames $K_1$ and $K_2$ in different directions relative to a third frame $K$ (not all coplanar) and to one another to place some constraints on it, via
$$
ds^2 = a(V_1)ds_1^2 = a(V_1) a(V_{12})ds_2^2 = a(V_2)ds_2^2 \ \to \ \ a(V_1) a(V_{12}) = a(V_2) \ \to \ \ a(V_{12}) = \frac{a(V_2)}{a(V_1)}.$$
Since the left-hand side $a(V_{12})$ depends on the angle between the velocities in $K_1$ and $K_2$, while the right-hand side doesn't, we see $a(V)$ must be a constant regardless of the argument $V$ thus $a^2 = a$ implies $a = 1$.
Now let's go to GR. Above we have discovered that the parameter $s'$ for which
$$ds_I'^2 = c^2 dt'^2 - dx'^2 - dy'^2 - dz'^2$$
holds in the frame $K'$ is such that
$$ds_I^2 = ds_I'^2$$
as long as $K$ and $K'$ are inertial frames.
Let's now compare to the case where we transform from an inertial frame to a non-inertial frame.
If we transform to a non-inertial frame, the usual first example being a frame $K'$ rotating about the z-axis with angular velocity $\Omega$ as
$$t = t', x = x' \cos (\Omega t) - y' \sin(\Omega t), y = x' \sin(\Omega t) + y' \cos(\Omega ), z = z'$$
then we have
\begin{align}
ds_I^2 &= c^2 dt^2 - dx^2 - dy^2 - dz^2 \\
&= [c^2 - \Omega^2(x'^2 + y'^2)]dt'^2 - dx'^2 - dy'^2 - dz'^2 + 2 \Omega y' dx' dt' - 2 \Omega x' dy' dt' \\
&= ds_I'^2 - \Omega^2(x'^2 + y'^2)dt'^2 + 2 \Omega y' dx' dt' - 2 \Omega x' dy' dt' \\
&= [ 1 - \Omega^2(x'^2 + y'^2)(\frac{dt}{ds_I'})'^2 + 2 \Omega y' (\frac{dx'}{ds_I'}) (\frac{dt'}{ds'}) - 2 \Omega x' \frac{dy'}{ds'} \frac{dt'}{ds_I'} ] ds_I'^2 \\
&= a(\Omega,t',x',...) ds_I'^2.
\end{align}
We now see that if I try to define the line element in a non-inertial frame in the exact same way as in the inertial case, I will find non-trivial $a = a(\Omega,t',x',...)$ factors. Since we are in a non-inertial frame, I no longer have symmetry arguments to make any general arguments about it's structure - I can't invoke homogeneity of space and time, nor isotropy of space, since these are not true in a non-inertial frame, so I can't say anything general about this coefficient, it depends on the properties of the non-inertial frame that I transform to. It's not not even true that $ds_I'^2$ needs to vanish for a light beam in $K'$ anymore, all I need are that $a(\Omega,...)ds_I'^2$ vanishes whenever $ds_I^2$ vanishes.
I just have to accept that the line element in a non-inertial frame takes a non-diagonal form where the coefficients can depend non-trivially on the coordinates, velocities etc... thus we define the line element in the non-inertial frame in this example as
$$
ds'^2 = [c^2 - \Omega^2(x'^2 + y'^2)]dt'^2 - dx'^2 - dy'^2 - dz'^2 + 2 \Omega y' dx' dt' - 2 \Omega x' dy' dt' = g_{\mu \nu}(x') dx'^{\mu} dx'^{\nu}.$$
When we study transformations to non-inertial reference frames, we usually consider it part of the study of general relativity and study transformations from an inertial to a non-inertial frame through the equality
$$ds_I^2 = \eta_{\mu \nu} dx^{\mu} dx^{\nu} = \eta_{\mu \nu} \frac{\partial x^{\mu}}{\partial x'^{\rho}} \frac{\partial x^{\nu}}{\partial x'^{\sigma}} dx'^{\rho} dx'^{\sigma} = g_{\rho \sigma}'(x') dx'^{\rho} dx'^{\sigma} = ds'^2$$
So in other words, I could write the above discussion as
$$ds_I^2 = (\frac{ds_I}{ds'})^2 ds'^2 = ds'^2 = (\frac{ds'}{ds_I'})^2 ds_I'^2.$$
It's clear that how I choose my parameter matters, and it's even more clear that in the non-inertial $K'$ frame that $ds_I'^2$ is different to $ds'^2$, again noting that even $ds'^2$ is also just a (scalar) parameter.
It is now hopefully clear that the closest Euclidean 3D analog to the above example would be something like setting $ds_I^2 = dx^2 + dy^2 + dz^2$ and $ds_I'^2 = dx'^2 + dy'^2 + dz'^2$ in an example with cross terms, e.g.
$$ds_I = \sqrt{dx^2 + dy^2 + dz^2} = \sqrt{dx'^2 + dy'^2 + dz'^2 + ...} = ds_I' \sqrt{1 + (.../ds_I'^2)} = ds'$$
and in general I can take any example, like the general one involving $g_{ij}$'s at the start of the post, and write $ds_I = \frac{ds_I}{ds_I'} ds_I' = ds' = \frac{ds'}{ds_I'} ds_I'$ as before.
All this is just the set-up to begin to understand the equivalence principle, that a particle moving in a gravitational field in an inertial frame is indistinguishable locally from a free particle moving in a non-inertial frame, thus a gravitational field is such that it affects the expression for $ds^2$ locally so that it takes on the form $ds^2 = g_{\mu \nu}(x) dx^{\mu} dx^{\nu}$. Globally however, I can simply invert the above transformation from a non-inertial frame and reclaim the Minkowski line element everywhere, but for an actual gravitational field this is not possible globally, indeed a general gravitational field depends on 10 independent functions, $g_{\mu \nu}$, a transformation to and from a non-inertial frame means only four of the $g_{\mu \nu}$ are independent.
For this argument please also see for example (among others)
