Closest points between two lines

The idea is that for the line segment of the shortest length, it has to be perpendicular to both the other lines.

Let the perpendicular line start from a point $P_1+t_1V_1$ of the first line, and have tangent vector $V_3$, i.e.: $$L_3=P_1+t_1V_1+t_3V_3$$

It should be that $$V_3\cdot V_2=0$$ $$V_3 \cdot V_1=0$$ which means that you can get $V_{3}$ as $$ V_3=V_2\times V_1$$

Now for this line to meet also the second line you need to have

$$P_1+t_1V_1+t_3V_3=P_2+t_2V_2$$

With this you have 3 linear equations in 3 variables, $t_1$,$t_2$ and $t_3$.

Once you solve for them, then:
- the distance between line 1 and line 2 will be $d=t_{3}/|V_{3}|$, to be taken, in case, as absolute value ($d$ is the "algebraic" distance in the direction of vector $V_{3}$);
- the closest points will of course be $Q_{1}=P_{1}+t_{1}V_{1}$ on line 1, and $Q_{2}=P_{2}+t_{2}V_{2}$ on line 2.

In case anyone else was, like me, looking to implement this, here is some example code in Python that follows @user2255770's answer:

from numpy import array, cross
from numpy.linalg import solve, norm

# define lines A and B by two points
XA0 = array([1, 0, 0])
XA1 = array([1, 1, 1])
XB0 = array([0, 0, 0])
XB1 = array([0, 0, 1])

# compute unit vectors of directions of lines A and B
UA = (XA1 - XA0) / norm(XA1 - XA0)
UB = (XB1 - XB0) / norm(XB1 - XB0)
# find unit direction vector for line C, which is perpendicular to lines A and B
UC = cross(UB, UA); UC /= norm(UC)

# solve the system derived in user2255770's answer from StackExchange: https://math.stackexchange.com/q/1993990
RHS = XB0 - XA0
LHS = array([UA, -UB, UC]).T
print(solve(LHS, RHS))
# prints "[ 0. -0.  1.]"