Post History

Perhaps the simplest way to get the formula is to think geometrically.

Let's say we wanted to reflect a point $x=(r\cos\theta,r\sin\theta)$ across the $x$-axis. In that case, we can simply negate $\theta$ giving $(r\cos(-\theta),r\sin(-\theta))=(r\cos\theta,-r\sin\theta)$ as expected.

If we want to reflect across a line at $\phi$ radians from the $x$-axis, we simply counter-rotate by $\phi$, i.e. subtract $\phi$ from the angles, to get into the above situation, negate to reflect, then rotate by $\phi$, i.e. add $\phi$ to the angles, to take things back to the original arrangement. Focusing just on the angles, this gives $-(\theta-\phi) + \phi = 2\phi - \theta$. The $\pi$ in the formulas in the OP is due to the fact that the $e'$ vector is the *normal* to the (hyper-)plane of reflection so $\phi$ is 90°; or $\frac{\pi}{2}$ radians off from it, i.e. $\phi = \theta_0 \pm \frac{\pi}{2}$ leading to $2\phi-\theta = 2\theta_0 \pm \pi - \theta$. (The sign of $\pi$ doesn't matter here.)

----

If you'd prefer an algebraic approach, here's an approach using [geometric algebra](https://web.archive.org/web/20230422191955/http://geocalc.clas.asu.edu/pdf/PrimerGeometricAlgebra.pdf).

First, we can represent the reflection of the vector $\mathbf x$ in the (hyper-)plane orthogonal to a unit vector $\mathbf e$ by the simple expression $-\mathbf{exe}$. Using the basic geometric algebra identities $\mathbf{uv} = \mathbf u\cdot\mathbf v + \mathbf u\wedge \mathbf v$ and $\mathbf u\cdot(\mathbf v\wedge \mathbf w) = (\mathbf u\cdot \mathbf v)\mathbf w - (\mathbf u\cdot \mathbf w)\mathbf v$ (and, implicitly, $\mathbf u\wedge \mathbf u = 0$), we get: $$\begin{align}
-\mathbf{exe}
&= -(\mathbf x\cdot \mathbf e)\mathbf e - \mathbf e(\mathbf x\wedge \mathbf e) \\\\
&= -(\mathbf x\cdot \mathbf e)\mathbf e - \mathbf e\cdot(\mathbf x\wedge \mathbf e) \\\\
&= -(\mathbf x\cdot \mathbf e)\mathbf e - (\mathbf e\cdot \mathbf x)\mathbf e + (\mathbf e\cdot \mathbf e)\mathbf x \\\\
&= \mathbf x - 2(\mathbf x\cdot \mathbf e)\mathbf e
\end{align}$$
reproducing the formula for $\sigma_{e'}$.

Now we need to talk about rotation. The general formula for rotating a vector $\mathbf x$ by a $\theta$ radians in the plane spanned by the unit vectors $\mathbf u$ and $\mathbf v$ is $e^{-B\theta/2}\mathbf xe^{B\theta/2}$ where $B=(\mathbf u\wedge \mathbf v)/\vert \mathbf u\wedge \mathbf v\vert$ and, since $B^2 = -1$, the exponential expands a la Euler's formula $e^{B\theta} = \cos\theta + B\sin\theta$. In the special case that $\mathbf x$ is in the span of $\mathbf u$ and $\mathbf v$, then $\mathbf x(\mathbf u\wedge \mathbf v) = -(\mathbf u\wedge \mathbf v)\mathbf x$ and the rotation formula simplifies to $e^{-B\theta}\mathbf x$. But if we compare this to $\mathbf{vu}$ we get: $$\mathbf{vu} = \mathbf v\cdot \mathbf u + \mathbf v\wedge \mathbf u = \mathbf v\cdot \mathbf u - \vert \mathbf u\wedge \mathbf v\vert B = \cos\theta - B\sin\theta = e^{-B\theta}$$ where $\theta$ is the angle between $\mathbf u$ and $\mathbf v$. In other words, $\mathbf{vu}$ represents the rotation that would rotate $\mathbf u$ into $\mathbf v$, evident from the fact that $\mathbf{vuu}=\mathbf v$.

Thus we can view $-\mathbf{exe}$ as "rotate $\mathbf e$ by the angle between $\mathbf e$ and $\mathbf x$ and then negate". In terms of angles, this means add the angle between $\mathbf e$ and $\mathbf x$, i.e. $\theta_0-\theta$, to the angle $\mathbf e$ is at, i.e $\theta_0$, and then negate which is equivalent to adding $\pm\pi$ to the angle. In symbols, $\theta_0 + (\theta_0-\theta) \pm \pi = 2\theta_0 - \theta \pm \pi$.

But screw intuition, I want to compute. Write $\mathbf x=e^{-B\theta}\mathbf{e}_0$ where $\mathbf{e}_0$ is the unit vector in the $x_0$ direction, and $\mathbf e = e^{-B\theta_0}\mathbf{e}_0$ where $B=\mathbf{e}_0\mathbf{e}_1$. We now get the above paragraph in symbols:
$$\begin{align}
-\mathbf{exe}
&= -e^{-B\theta_0}\mathbf{e}_0 e^{-B\theta}\mathbf{e}_0 e^{-B\theta_0}\mathbf{e}_0 \\\\
&= -e^{-B\theta_0}e^{B\theta}\mathbf{e}_0\mathbf{e}_0 e^{-B\theta_0}\mathbf{e}_0 \tag{$\mathbf{e}_0^2 = 1$}\\\\
&= -e^{-B(\theta_0-\theta)}e^{-B\theta_0}\mathbf{e}_0 \\\\
&= -e^{-B(2\theta_0-\theta)}\mathbf{e}_0 \\\\
&= e^{\mp B\pi}e^{-B(2\theta_0-\theta)}\mathbf{e}_0 \tag{$-1 = e^{\mp B\pi}$} \\\\
&= e^{-B(2\theta_0-\theta\pm\pi)}\mathbf{e}_0 \\\\
&= \cos(2\theta_0-\theta\pm\pi)\mathbf{e}_0 + \sin(2\theta_0-\theta\pm\pi)\mathbf{e}_1
\end{align}$$

----

As for matrices, you can recover them if you want by introducing a basis and seeing how each basis vector is transformed. But there's no reason to do this. Talking about matrices (as opposed to linear transformations) necessarily implies introducing a basis and coordinates, and that's simply not necessary most of the time. Often it will be necessary to introduce a basis to concretely specify a vector, but the formulas can be manipulated and derived without doing that. Usually the formulas are simpler and more general by avoiding coordinates. You can see above that the only time I introduced a basis in the second part was to connect to the coordinate-based expressions in the question.