Of course this theorem has a geometric interpretation! In a sense, it's a multidimensional analogue of «the volume of a parallelepiped is the product of the area of its base and its height».
3. Let's start with $3\times3$ case:
$$
\left|\begin{matrix}u_1&u_2&u_3\\v_1&v_2&v_3\\w_1&w_2&w_3\end{matrix}\right|=
u_1\left|\begin{matrix}v_2&v_3\\w_2&w_3\end{matrix}\right|
-u_2\left|\begin{matrix}v_1&v_3\\w_1&w_3\end{matrix}\right|
+u_3\left|\begin{matrix}v_1&v_2\\w_1&w_2\end{matrix}\right|.
$$
LHS is the volume of the parallelepiped spanned by three vectors, $u$, $v$ and $w$. What's the meaning of RHS? Clearly that's a scalar product of $u$ with something — namely, with the vector
$$
\left(\left|\begin{matrix}v_2&v_3\\w_2&w_3\end{matrix}\right|,
-\left|\begin{matrix}v_1&v_3\\w_1&w_3\end{matrix}\right|,\left|\begin{matrix}v_1&v_2\\w_1&w_2\end{matrix}\right|\right)=
\left|\begin{matrix}\overrightarrow{e_1}&\overrightarrow{e_2}&\overrightarrow{e_3}\\v_1&v_2&v_3\\w_1&w_2&w_3\end{matrix}\right|
$$ — i.e. with vector product of $v$ and $w$.
So the formula we get is $vol\langle u,v,w\rangle=(u,[v,w])$; now by the (geometrical) definition of scalar product it's $area\langle v,w\rangle\cdot (|u|\cdot\sin\phi)$, and the first factor is the area of the base and the second one is the height of our parallelepiped.
n. Consider the (general) case of vectors in $n$-dimensional space $V$. In RHS of the theorem we again see a scalar product of the first vector, $v$, with a vector $B$ (in coordinate-free language it really lives in $\Lambda^{n-1}V$, but let's ignore this for now) with coordinates $C_{1i}$.
The question is, what is the geometric meaning of $B$. Let me give 3 (closely related) answers.
- By the very same cofactor theorem it measures the [(n-1)-dimensional] area of projection of the base of our $n$-parallelepiped (i.e. $(n-1)$-parallelepiped spanned by all vectors but $v$) on different hyperplanes; more precisely, the area of the projection on the hyperplanes orthogonal to a unit vector $v$ is the scalar product $(B,v)$.
- Let's prove the cofactor theorem instead of using it. The function $(B,x)$ is linear in $x$. For a basis vector $x=e_i$ we have $(B,x)=C_{1i}$, which (up to sign, at least) is the area of the span of projections of our vectors on the hyperplane orthogonal to $e_i$. So $(B,x)$ is indeed the area of the projection of the base on the hyperplane orthogonal to $x$ (multiplied by $|x|$ and taken with appropriate signs).
- Even better, since everything is invariant under (special) orthogonal transforms, let's change basis to make $v$ a scalar multiple of $e_1$. Now the statement «$(B,v)$ is the $|v|$ times the area of the projection» became obvious (we literally multiply $|v|$ by the cofactor manifestly equal to this area — well, it was discussed in (2) anyway).
Now I must admit the statement we get is more like «the volume of a parallelepiped $\langle u,\text{base}\rangle$ is the product of the length of $u$ and the area of the projection of its base on the hyperplane orthogonal to $u$» — but it's of course equivalent to «the volume of a parallelepiped is the product of the area of its base and its height».
I believe the recurrence should be $f(n)=n(f(n-1)+1)$ with $f(2)=2.$ You can absorb the $(-1)^{j+k}$ in the addition of the products, so each term just takes one multiply. This gives $0,2,9,40,205,1236$, which is OEIS A038156 and matches your author's expression except that the upper limit of the sum should be $n-1$, not $n$. Another formula is $f(n)=\lfloor n!(e-1) \rfloor -1$
We can verify your author's formula from the recurrence by induction. The base case $f(2)=2=2!\frac 1{1!}$ works. Now assume we have $f(n)=n!\sum_{k=1}^{n-1}\frac 1{k!}$ We have $$\begin {align} f(n+1) &=(n+1)\left(f(n)+1\right)\\
&=(n+1)\left(n!\sum_{k=1}^{n-1}\frac 1{k!}+1\right)\\
&=(n+1)!\left(\sum_{k=1}^{n-1}\frac 1{k!}+ \frac 1{n!}\right)\\
&==(n+1)!\left(\sum_{k=1}^{n}\frac 1{k!}\right) \end {align}$$
Best Answer
Here is a non-clever approach that is still better than brute force. Subtract the bottom row from each of the preceding rows. This does not change the determinant. The result is $$\begin{pmatrix} 0 & x-w & x^2-w^2 & x^3-w^3 \\ 0 & y-w & y^2-w^2 & y^3-w^3 \\ 0 & z-w & z^2-w^2 & z^3-w^3 \\ 1 & w & w^2 & w^3 \end{pmatrix}$$ Expanding along the first column, find that the determinant is $$-\det \begin{pmatrix} x-w & x^2-w^2 & x^3-w^3 \\ y-w & y^2-w^2 & y^3-w^3 \\ z-w & z^2-w^2 & z^3-w^3 \end{pmatrix}$$ Factor out $(x-w)(y-w)(z-w)$: $$-(x-w)(y-w)(z-w) \det \begin{pmatrix} 1 & x+w & x^2+xw+w^2 \\ 1 & y+w & y^2+yw+w^2 \\ 1 & z+w & z^2+zw+w^2 \end{pmatrix}$$ Again subtract the bottom row from the first two. $$-(x-w)(y-w)(z-w) \det \begin{pmatrix} 0 & x-z & x^2-z^2+(x-z)w \\ 0 & y-z & y^2-z^2+(y-z)w \\ 1 & z+w & z^2+zw+w^2 \end{pmatrix}$$ This becomes a $2\times 2$ determinant, from which you can factor $(x-z)(y-z)$: $$-(x-w)(y-w)(z-w)(x-z)(y-z) \det \begin{pmatrix} 1 & x+z+1 \\ 1 & y+z+1 \end{pmatrix}$$ and conclude.