Another reference which takes a very explicit point of view (i.e. entirely matrix calculations and explicit circle inversions) is Beardon's The Geometry of Discrete Groups, which deduces the invariance in $n$ dimensions from the following useful formula: if $\sigma$ is the reflection in the Euclidean sphere of radius $r$ and centre $a$, then

$$ |\sigma(y) - \sigma(x)| = \frac{r^2 |y-x|}{|x-a||y-a|} $$

(this is formula 3.1.5 on p.26). Observe now that the extension $\tilde{\phi}$ to $\mathbb{R}^{n+1} $ of any inversion $\phi$ in $ \mathbb{R}^n$ is a reflection in some circle of radius $r$ centred on $ a $ with $a_{n+1} = 0 $, and explicitly writing down the formula for circle inversion gives that the $(n+1)$th coordinate of $\tilde{\phi}(x) $ is

$$ \frac{r^2 x_{n+1}}{|x-a|^2}. $$

But this formula is almost $ x_{n+1} $ times the first displayed formula above, and this similarity suggests

\begin{align*} \frac{|\tilde{\phi}(y)-\tilde{\phi}(x)|^2}{\tilde{\phi}(x)_{n+1} \tilde{\phi}(y)_{n+1}} &= \frac{r^4 |y-x|^2}{|x-a|^2|y-a|^2 \frac{r^2 x_{n+1}}{|x-a|^2} \frac{r^2 y_{n+1}}{|y-a|^2}}\\ &= \frac{|y-x|^2}{x_{n+1} y_{n+1}} \end{align*}

i.e. the form $ |y-x|^2/x_{n+1} y_{n+1} $ is invariant under Poincare extensions of sphere inversions. The same is true for extensions of plane reflections, and since Poincare extension is a homomorphism we get that the extension of every Moebius transformation preserves that form; but the form is just the hyperbolic metric:

$$ d(x,y) = \int_\gamma \frac{|y-x|^2}{x_{n+1} y_{n+1}} ds. $$

This calculation is found on pp.34--35 of Beardon (I have added a very small amount of detail). At this point of the book he has only done the structure of Moebius maps as products of reflections, i.e. it is very hands-on and requires essentially no "heavy" machinery. The downside of this approach is that the result appears magically out of thin air as opposed to other appproaches which come from symmetric spaces and Lie groups, e.g. showing that $ \mathrm{PSL}(2,\mathbb{C}) $ is naturally isomorphic to $ O(2,1) $ as a Lie group and hence admits an action by isometries on $ \mathbb{H}^3 $ which you then show is the usual action by fractional linear transformations.

Another reference which takes a very explicit point of view (i.e. entirely matrix calculations and explicit circle inversions) is Beardon's The Geometry of Discrete Groups, which deduces the invariance in $n$ dimensions from the following useful formula: if $\sigma$ is the reflection in the Euclidean sphere of radius $r$ and centre $a$, then we can estimate the derivative of $ \sigma $ by

$$ \frac{|\sigma(y) - \sigma(x)|}{|y-x|} = \frac{r^2}{|x-a||y-a|} $$

(this is formula 3.1.5 on p.26). Observe now that the extension $\tilde{\phi}$ to $\mathbb{R}^{n+1} $ of any inversion $\phi$ in $ \mathbb{R}^n$ is a reflection in some circle of radius $r$ centred on $ a $ with $a_{n+1} = 0 $, and explicitly writing down the formula for circle inversion gives that the $(n+1)$th coordinate of $\tilde{\phi}(x) $ is

$$ \frac{r^2 x_{n+1}}{|x-a|^2}. $$

But this formula is almost $ x_{n+1} $ times the first displayed formula above, and this similarity suggests

\begin{align*} \frac{|\tilde{\phi}(y)-\tilde{\phi}(x)|^2}{\tilde{\phi}(x)_{n+1} \tilde{\phi}(y)_{n+1}} &= \frac{r^4 |y-x|^2}{|x-a|^2|y-a|^2 \frac{r^2 x_{n+1}}{|x-a|^2} \frac{r^2 y_{n+1}}{|y-a|^2}}\\ &= \frac{|y-x|^2}{x_{n+1} y_{n+1}} \end{align*}

i.e. the form $ |y-x|^2/x_{n+1} y_{n+1} $ is invariant under Poincare extensions of sphere inversions. The same is true for extensions of plane reflections, and since Poincare extension is a homomorphism we get that the extension of every Moebius transformation preserves that form; but the form is just the hyperbolic metric:

$$ \sum_{k=0}^{N} \sqrt{\frac{|\gamma(t^k) - \gamma(t^{k-1})|^2}{\gamma(t^k)_{n+1} \gamma(t^{k-1})_{n+1}}} \to \int \frac{|d\gamma(t)|}{\gamma(t)_{n+1}} $$

(where I have written an approximation to the length of the curve $ \gamma $ by taking a partition $ t^0, \dots, t^N $: the function under the sum is clearly the square root of the invariant form above, and it approaches the integral at the right where the integrand is the usual expression for the hyperbolic metric $|dx|/x_{n+1} $).

This calculation is found on pp.34--35 of Beardon (I have added a very small amount of detail). At this point of the book he has only done the structure of Moebius maps as products of reflections, i.e. it is very hands-on and requires essentially no "heavy" machinery. The downside of this approach is that the result appears magically out of thin air as opposed to other appproaches which come from symmetric spaces and Lie groups, e.g. showing that $ \mathrm{PSL}(2,\mathbb{C}) $ is naturally isomorphic to $ O(2,1) $ as a Lie group and hence admits an action by isometries on $ \mathbb{H}^3 $ which you then show is the usual action by fractional linear transformations.

Another reference which takes a very explicit point of view (i.e. entirely matrix calculations and explicit circle inversions) is Beardon's The Geometry of Discrete Groups, which deduces the invariance in $n$ dimensions from the following useful formula: if $\sigma$ is the reflection in the Euclidean sphere of radius $r$ and centre $a$, then

$$ |\sigma(y) - \sigma(x)| = \frac{r^2 |y-x|}{|x-a||y-a|} $$

(this is formula 3.1.5 on p.26). Observe now that the extension $\tilde{\phi}$ to $\mathbb{R}^{n+1} $ of any inversion $\phi$ in $ \mathbb{R}^n$ is a reflection in some circle of radius $r$ centred on $ a $ with $a_{n+1} = 0 $, and explicitly writing down the formula for circle inversion gives that the $(n+1)$th coordinate of $\tilde{\phi}(x) $ is

$$ \frac{r^2 x_{n+1}}{|x-a|^2}. $$

But this formula is just $ x_{n+1} $ times the first displayed formula above, and this similarity suggests

\begin{align*} \frac{|\tilde{\phi}(y)-\tilde{\phi}(x)|^2}{\tilde{\phi}(x)_{n+1} \tilde{\phi}(y)_{n+1}} &= \frac{r^4 |y-x|^2}{|x-a|^2|y-a|^2 \frac{r^2 x_{n+1}}{|x-a|^2} \frac{r^2 y_{n+1}}{|y-a|^2}}\\ &= \frac{|y-x|^2}{x_{n+1} y_{n+1}} \end{align*}

i.e. the form $ |y-x|^2/x_{n+1} y_{n+1} $ is invariant under Poincare extensions of sphere inversions. The same is true for extensions of plane reflections, and since Poincare extension is a homomorphism we get that the extension of every Moebius transformation preserves that form; but the form is just the hyperbolic metric:

$$ d(x,y) = \int_\gamma \frac{|y-x|^2}{x_{n+1} y_{n+1}} ds. $$

This calculation is found on pp.34--35 of Beardon (I have added a very small amount of detail). At this point of the book he has only done the structure of Moebius maps as products of reflections, i.e. it is very hands-on and requires essentially no "heavy" machinery. The downside of this approach is that the result appears magically out of thin air as opposed to other appproaches which come from symmetric spaces and Lie groups, e.g. showing that $ \mathrm{PSL}(2,\mathbb{C}) $ is naturally isomorphic to $ O(2,1) $ as a Lie group and hence admits an action by isometries on $ \mathbb{H}^3 $ which you then show is the usual action by fractional linear transformations.

Another reference which takes a very explicit point of view (i.e. entirely matrix calculations and explicit circle inversions) is Beardon's The Geometry of Discrete Groups, which deduces the invariance in $n$ dimensions from the following useful formula: if $\sigma$ is the reflection in the Euclidean sphere of radius $r$ and centre $a$, then

$$ |\sigma(y) - \sigma(x)| = \frac{r^2 |y-x|}{|x-a||y-a|} $$

(this is formula 3.1.5 on p.26). Observe now that the extension $\tilde{\phi}$ to $\mathbb{R}^{n+1} $ of any inversion $\phi$ in $ \mathbb{R}^n$ is a reflection in some circle of radius $r$ centred on $ a $ with $a_{n+1} = 0 $, and explicitly writing down the formula for circle inversion gives that the $(n+1)$th coordinate of $\tilde{\phi}(x) $ is

$$ \frac{r^2 x_{n+1}}{|x-a|^2}. $$

But this formula is almost $ x_{n+1} $ times the first displayed formula above, and this similarity suggests

\begin{align*} \frac{|\tilde{\phi}(y)-\tilde{\phi}(x)|^2}{\tilde{\phi}(x)_{n+1} \tilde{\phi}(y)_{n+1}} &= \frac{r^4 |y-x|^2}{|x-a|^2|y-a|^2 \frac{r^2 x_{n+1}}{|x-a|^2} \frac{r^2 y_{n+1}}{|y-a|^2}}\\ &= \frac{|y-x|^2}{x_{n+1} y_{n+1}} \end{align*}

i.e. the form $ |y-x|^2/x_{n+1} y_{n+1} $ is invariant under Poincare extensions of sphere inversions. The same is true for extensions of plane reflections, and since Poincare extension is a homomorphism we get that the extension of every Moebius transformation preserves that form; but the form is just the hyperbolic metric:

$$ d(x,y) = \int_\gamma \frac{|y-x|^2}{x_{n+1} y_{n+1}} ds. $$

This calculation is found on pp.34--35 of Beardon (I have added a very small amount of detail). At this point of the book he has only done the structure of Moebius maps as products of reflections, i.e. it is very hands-on and requires essentially no "heavy" machinery. The downside of this approach is that the result appears magically out of thin air as opposed to other appproaches which come from symmetric spaces and Lie groups, e.g. showing that $ \mathrm{PSL}(2,\mathbb{C}) $ is naturally isomorphic to $ O(2,1) $ as a Lie group and hence admits an action by isometries on $ \mathbb{H}^3 $ which you then show is the usual action by fractional linear transformations.