Why is the scalar product of two four-vectors Lorentz-invariant?

Frankly, you're looking at this backwards.

Why is $A^\mu g_{\mu\nu}B^{\nu}$ Lorentz-invariant?

That's the wrong way around: $A^\mu g_{\mu\nu}B^{\nu}$ is Lorentz invariant because Lorentz transformations are defined as the class of transformations that leaves $A^\mu g_{\mu\nu}B^{\nu}$ invariant.

Generally, if you transform $A^\mu$ and $B^\mu$ by some linear transformation with the transformation matrix $\Lambda^\mu_{\ \ \nu}$, then their transformed values will be $\tilde A^\mu = \Lambda^\mu_{\ \ \nu}A^\nu$ and $\tilde B^\mu = \Lambda^\mu_{\ \ \nu} B^\nu$ (using Einstein summations). This means that the transformed inner product will be \begin{align} \tilde A^\mu g_{\mu\nu} \tilde B^{\nu} & = (\Lambda^\mu_{\ \ \alpha}A^\alpha) g_{\mu\nu} (\Lambda^\nu_{\ \ \beta}B^\beta) \\ & = A^\alpha (\Lambda^\mu_{\ \ \alpha} g_{\mu\nu} \Lambda^\nu_{\ \ \beta} )B^\beta \\ & = A^\mu (\Lambda^\gamma_{\ \ \mu} g_{\gamma\delta} \Lambda^\delta_{\ \ \nu} )B^\nu \qquad\qquad\text{(by re-labelling)} \\ & \stackrel{\text{required}}= A^\mu g_{\mu\nu} B^\nu. \end{align} Thus, for $A^\mu g_{\mu\nu} B^\nu$ to be invariant we require that $$ A^\mu (\Lambda^\gamma_{\ \ \mu} g_{\gamma\delta} \Lambda^\delta_{\ \ \nu} )B^\nu = A^\mu g_{\mu\nu} B^\nu $$ for all $A^\mu$ and $B^\mu$, and by judicious choices of those vectors (basically running each independently over the basis in use) that can only be the case if $$ \Lambda^\gamma_{\ \ \mu} g_{\gamma\delta} \Lambda^\delta_{\ \ \nu} = g_{\mu\nu}, $$ which forms the core requirement on $\Lambda^\mu_{\ \ \nu}$ for it to be a Lorentz transformation.

Here's the way to think about this -- why is the standard Euclidean dot product, $\sum x_iy_i$ interesting? Well, it is interesting primarily from the perspective of rotations, due to the fact that rotations leave dot products invariant. The reason this is so is that this dot product can be written as $|x||y|\cos\Delta\theta$, and rotations leave magnitudes and relative angles invariant.

Is the standard Euclidean norm $|x|$ invariant under Lorentz transformations? Of course not -- for instance, $\Delta t^2+\Delta x^2$ is clearly not invariant, but $\Delta t^2-\Delta x^2$ is. Similarly, $E^2+p^2$ is not important, but $E^2-p^2$ is. The reason this is the case is that Lorentz boosts are fundamentally skew transformations, which means the invariant locus is a hyperbola, not a circle. So you have $\cosh^2 \xi - \sinh^2 \xi = 1$, and $x_0^2-x_1^2$ is the right way to think of the norm on Minkowski space.

Similarly, Lorentz boosts change the rapidity $\xi$ by a simple displacement, so $\Delta \xi$ is invariant. From this point, it's a simple exercise to show that

$$|x||y|\cosh\xi=x_0y_0-x_1y_1$$

(as for the remaining dimensions -- remember that the standard Euclidean dot product is still relevant in space, so you just need to write $x_0y_0-x\cdot y=x_0y_0-x_1y_1-x_2y_2-x_3y_3$.)