Derive Distance of skew lines with methods from analysis It is well known, that given two skew lines
\begin{align*}
g: x &= a + \lambda b \quad \\
h: y &= c + \mu d
\end{align*}
the distance is given by 
$$D = \frac{|(a - c) \cdot b \times d|}{|b \times d|}$$
See e.g. here
I thought, it would be a nice exercise to derive the distance using an optimization 
method, but I get stuck. Short version: Why is
\begin{equation*}
 \begin{aligned}
   \left\| \left(a - c\right) - \frac{|b|^2\cdot |d|^2}{|b \times d|^2} \left[\vphantom{\frac{b}{b}}\right.\right.
   & \left.\left. \left(\left(a - c\right)\cdot \left(\frac{b}{|b|} - \frac{b \cdot d}{|b|\cdot|d|} \cdot \frac{d}{|d|}\right)\right) \cdot \frac{b}{|b|} \right.\right. + \\
   & \left.\left. \left(\left(a - c\right)\cdot \left(\frac{d}{|d|} - \frac{b \cdot d}{|b|\cdot|d|} \cdot \frac{b}{|b|}\right)\right) \cdot \frac{d}{|d|} 
   \right] \right\| = \frac{|(a - c) \cdot b \times d|}{|b \times d|}
 \end{aligned}
\end{equation*}
It seems clear to me, that if you subtract the orthogonal projections of $a - c$ in
direction of $b$ and $d$, the result must be proportional to $b \times d$. 
On the other hand, there must be a clear line of arguments, why 
$\frac{|b|^2\cdot |d|^2}{|b \times d|^2}$ exactly does the job.

Long version. For notation and reference see e.g. here: 
The difference vector between two of those skew lines is
$$D = |y - x| = |c - a + \mu d - \lambda b|$$
and the local extrema $(\lambda, \mu)$ is found to be at
$$\frac{1}{\mbox{det } A}
\left(\begin{matrix}
 d^2 & d\cdot b \\
 d \cdot b & b^2  
\end{matrix}\right) \cdot \left(\begin{matrix}
 (c-a) \cdot b \\
 -(c-a) \cdot d  
\end{matrix}\right)$$
where $$\mbox{det A} = b^2 d^2 - (d \cdot b)^2 = |b \times d|^2$$
If you insert above solution into $D$, you get 
\begin{equation*}
 \begin{aligned}
   D = |c - a + \mu d - \lambda b|
 \end{aligned}
\end{equation*}
\begin{equation*}
 \begin{aligned}
   D = \left\| \left(a - c\right) - \left[\vphantom{\frac{b}{b}}\right.\right.
   & \frac{|d|^2}{|b \times d|^2} \left.\left. \left(\left(a - c\right)\cdot \left(b - \frac{b \cdot d}{|b|\cdot|d|^2} \cdot d\right)\right) \cdot b \right.\right.+\\
   & \frac{|b|^2}{|b \times d|^2} \left.\left. \left(\left(a - c\right)\cdot \left(d - \frac{b \cdot d}{|b|^2\cdot|d|} \cdot b\right)\right) \cdot d 
   \right] \right\|
 \end{aligned}
\end{equation*}
\begin{equation*}
 \begin{aligned}
   D = \left\| \left(a - c\right) - \frac{|b|^2\cdot |d|^2}{|b \times d|^2} \left[\vphantom{\frac{b}{b}}\right.\right.
   & \left.\left. \left(\left(a - c\right)\cdot \left(\frac{b}{|b|} - \frac{b \cdot d}{|b|\cdot|d|} \cdot \frac{d}{|d|}\right)\right) \cdot \frac{b}{b} \right.\right.+\\
   & \left.\left. \left(\left(a - c\right)\cdot \left(\frac{d}{|d|} - \frac{b \cdot d}{|b|\cdot|d|} \cdot \frac{b}{|b|}\right)\right) \cdot \frac{d}{d} 
   \right] \right\|
 \end{aligned}
\end{equation*}
Then I don't see how to proceed...
As pointed out by @ancientmathematician, we could set $u=a-c$, and assume $|b|=|d|=1$ without loss of generality, then
\begin{equation*}
 \begin{aligned}
   D = \left\| u - \frac{1}{|b \times d|^2} \left[\vphantom{\frac{b}{b}}\right.\right.
   & \left.\left. \left(u\cdot \left(b - (b \cdot d) \cdot d\right)\right) \cdot b \right.\right.+\\
   & \left.\left. \left(u\cdot \left(d - (b \cdot d) \cdot b\right)\right) \cdot d \vphantom{\frac{b}{b}} \right] \right\|
 \end{aligned}
\end{equation*}
and $|b \times d|^{-2}= \sin(\alpha)^{-2}$. But the occurence of the angle between $b$ and $d$ at that place is quite puzzling to me. Do I miss something obvious?
Initially I thought
\begin{equation*}
 \begin{aligned}
   D = \left\| u - \left[\right.\right.
   \left.\left. \left(u\cdot \left(b - (b \cdot d) \cdot d\right)\right) \cdot b \right.\right.+
   \left.\left. \left(u\cdot \left(d - (b \cdot d) \cdot b\right)\right) \cdot d 
   \right] \right\|
 \end{aligned}
\end{equation*}
which is wrong...

Without loss of generality, let $\ c=0,\ $ but we do not let $\ |b|=|d|=1\ $
 because we would lose the valuable homogeneous property. 
We recall the standard scalar
 Triple product identity:
$$ \Delta := (a \cdot (b \times d))^2 = \begin{vmatrix} a\cdot a & a \cdot b & a\cdot d \\
b \cdot a & b\cdot b & b\cdot d \\
d \cdot a & d\cdot b & d\cdot d 
 \end{vmatrix}. \tag{1} $$ 
$$ \textrm{Let }\quad b_v := \frac{b}{|b|} \!-\! \frac{b\cdot d}{|b|\,|d|} \frac{d}{|d|},
   \quad   d_v := \frac{d}{|d|} \!-\! \frac{b\cdot d}{|b|\,|d|} \frac{b}{|b|},
  \quad \textrm{ and} \tag{2} $$
$$ a_v := a\ |b\times d|^2 - |b|^2|d|^2 \Big(
 a\cdot b_v\frac{b}{|b|} + a\cdot d_v\frac{d}{|d|}\Big),
 \quad D := a_v/|b\times d|^2. \tag{3}  $$
A bit of simplification of equation $(2)$ gives us
$$ |b|\ |d|^2\ b_v = b\ |d|^2 - d\ (b\cdot d), \quad
   |b|^2\ |d|\ d_v = d\ |b|^2 - b\ (b\cdot d). \tag{4} $$
Substituting equation $(3)$ in equation $(4)$ gives us
$$ a_v = (a_0)a + (b_0)b + (d_0)d \quad \textrm{ with } \tag{5} $$
$$ a_0 :=|b\!\times\!d|^2,\quad 
b_0 := -a\!\cdot\!\big(b\ |d|^2 \!-\! d\ (b\!\cdot\! d)\big),\quad
d_0 := -a\!\cdot\!\big(d\ |b|^2 \!-\! b\ (b\!\cdot\! d)\big). \tag{6} $$
Taking the dot product of equation $(5)$ we get
$$ a_v\!\cdot\!a_v = (a\!\cdot\!a)a_0^2 \!+\! (b\!\cdot\!b)b_0^2 
  \!+\! (d\!\cdot\!d)d_0^2 \!+\! 2(a\!\cdot\!b)a_0b_0 \!+\!
   2(a\!\cdot\!d)a_0d_0 \!+\! 2(b\!\cdot\!d)b_0d_0. \tag{7}$$
Using $\ |b\!\times\!d|^2 = |b|^2|d|^2\!-\!(b\!\cdot\!d)^2,\ $
 and equations $(1)$ and $(6)$ and factoring equation $(7)$ we get
$$ |a_v|^2 = a_v\cdot a_v = |b\times d|^2 \Delta =
   |b\times d|^2 |a\cdot (b\times d)|^2 \tag{8}$$
Solving for $\ |a_v|\ $ and using equation $(3)$ we get
$$ |a_v| = |b\times d|\ |a\cdot b\times d|,\quad
   |D| = |a\cdot b\times d|/|b\times d| \tag{9} $$
which is what we wanted to prove.
