Solution (source code)

= Solution

Write $\mathbf w=\mathbf a-\mathbf c$, $B=\mathbf b\cdot\mathbf b$, $D=\mathbf d\cdot\mathbf d$ and $E=\mathbf b\cdot\mathbf d$. Differentiating $|\mathbf w+\lambda\mathbf b-\mu\mathbf d|^2$ with respect to the two parameters gives
$$
B\lambda-E\mu=-\mathbf w\cdot\mathbf b,\qquad
E\lambda-D\mu=-\mathbf w\cdot\mathbf d.
$$
For nonparallel directions, $BD-E^2=|\mathbf b\times\mathbf d|^2>0$. Solving this <linear system> yields
$$
\boxed{\lambda=\frac{(\mathbf b\cdot\mathbf d)[(\mathbf a-\mathbf c)\cdot\mathbf d]-|\mathbf d|^2[(\mathbf a-\mathbf c)\cdot\mathbf b]}{|\mathbf b|^2|\mathbf d|^2-(\mathbf b\cdot\mathbf d)^2}.}
$$
The quadratic being minimized is strictly convex in the two parameters under this nonparallel assumption, so the critical point gives the unique pair of closest points. Parallel lines do not select a unique value of $\lambda$.