For four non-parallel, non-intersecting lines in 3D space, two other lines intersect with all four lines simultaneously. How are they found? Using CAD software, I have been able to make a construction in which four arbitrary lines in space (infinitely long) have two lines that intersect with all four. There are only ever two solutions that the software finds and degrees of freedom analysis vaguely suggests this is the expected result.
So we have four lines in space:
$$ \vec I = \vec I_0 + \alpha \vec I_1$$
$$\vec J = \vec J_0 + \beta \vec J_1$$
$$ \vec K = \vec K_0 + \gamma \vec K_1$$
$$ \vec L = \vec L_0 + \kappa \vec L_1$$
where $\vec \cdot$ denotes a vector. Two more lines make the intersections:
$$ \vec X = \vec X + \eta \vec X_1$$
$$ \vec Y = \vec Y_0 + \zeta \vec Y_1$$
There will be a unique set of $\alpha$, $\beta$, $\gamma$, $\kappa$ and $\eta$ where the line $X = I$, $X = J$, $X = K$ and $X = L$
There will be another unique set of $\alpha$, $\beta$, $\gamma$, $\kappa$ and $\zeta$ where the line $Y = I$, $Y = J$, $Y = K$ and $Y = L$.
I am supposing the six scalars will each have a parabolic function of a common parameter, say $f(t)$, where the solutions to $f(t) = 0$ are the unique values.
This also creates a challenge of identifying which solution corresponds to which intersection.
I can vaguely see a way through the mathematics but the complexity is getting beyond me. Can anyone help?
 A: Consider any three of the lines, say $I,J,K$. These three lines lie on a ruling of a unique quadric $Q$, which is either a hyperboloid of one sheet in the generic case, or a hyperbolic paraboloid if all three lines are parallel to the same plane. Any line that intersects all three lines is part of the opposite ruling of $Q$. Conversely, any point on $Q$ determines a line intersecting $I,J,K$, namely the member of the opposite ruling passing through said point.

Therefore, the lines intersecting all four lines $I,J,K,L$ are determined by the points of intersection of $Q$ and $L$. Each point of intersection $p$ gives a line which intersects $L$ (because they both pass through $p$) and intersects $I,J,K$ (because it is a member of the opposite ruling of $Q$).
Generically, there are either zero or two intersections. It is possible to have one intersection if $L$ is tangent to $Q$, but also in other cases -- for example, if $L$ is parallel to an asymptote(?) of the hyperboloid (or is that also a tangent to $Q$, in a projective-geometry sense?).
A: Here is an explanation of why there should be two solutions in the generic case.
First consider the problem for only three lines.
If you take three lines whose directions are not coplanar, who don't already intersect each other, and with no parallel pairs, then you can move and rescale the $x,y,z$-axis such that the equations of the lines are given by $L_1 : (x=0,y=1), L_2 : (y=0,z=1), L_3 : (x=1,z=0)$.
Three points $(x_i,y_i,z_i)$ are colinear if and only if the space spanned by the three $4$-vectors $(x_i,y_i,z_i,1)$ has rank $ \le 2$.
So, in our case, a point $(0,1,z) \in L_1$ forms a line with points $(x,0,1) \in L_2$ and $(1,y,0) \in L_3$ when
$rank\begin{pmatrix}x & 0 & 1 & 1 \\ 1 & y & 0 & 1 \\ 0 & 1 & z & 1 \end{pmatrix} \le 2$
Looking at the last three columns, the determinant has to vanish, so that  
$0 = \det\begin{pmatrix} 0 & 1 & 1 \\ y & 0 & 1 \\ 1 & z & 1 \end{pmatrix} = 1+yz-y$
And so $y = 1/(1-z)$. And for symmetry reasons, we also have $z-zx-1 = 0$ so that $x=1-1/z$.
Therefore, given one point $(0,1,z)$ on $L_1$, if $z \neq 0,1$ then there is a unique line going through it and that also intersects $L_2$ and $L_3$ (And you can let go of that restriction if you work in projective space and projective lines).
And $y$ (and $x$) are obtained from $z$ by applying a degree $1$ rational fraction.
If you are a point sitting on $L_1$ and you are looking around twisting your head in all directions, you "see" $L_2$ and $L_3$ "intersect" from your point of view at only one point.

Taking the reunion of this one-dimensional family of lines gives a surface (a ruled quadratic surface), and if you have a fourth line $L_4$, determining the two intersecting lines is the same thing as determining the intersection of $L_4$ with the surface.
To do this you want an equation of the surface. Given a point $(0,1,t)$ on $L_1$, you get the line going through $(0,1,t),(1- \frac 1t,0,1)$ and $(1, \frac 1{1-t},0)$. Choosing a parametrisation of those lines gives you a parametrisation of the surface, for example
$$s(1,\frac 1{1-t},0)+(1-s)(0,1,t) = (s,1-s+ \frac s{1-t},(1-s)t)$$
Then you have to work a bit to get the equation :
if $(x,y,z) = (s,1-s+\frac s{1-t},(1-s)t)$ then
$s = x$
$t = \frac z {1-s} = \frac z{1-x}$
so that :
$y = 
1-x+\frac x{1-\frac z{1-x}} = 1-x+\frac {x(1-x)}{1-x-z} = \frac{1-x-z+xz}{1-x-z}$.
From there you get the equation
$y(1-x-z) = 1-x-z+xz$
or $x+y+z = 1+xy+yz+xz$.
As promised this is a degree $2$ equation so it intersects most lines at two points.

Alternatively, forget about how we picked a very special reference frame. Given any affine parametrization of $L_1$ and $L_2$ (your $I = I_0 + \alpha I_1$ for points on $L_1$ and $J = J_0 + \beta J_1$ for points on $L_2$) there will be a degree $1$ rational fraction $f$ such that $(IJ)$ intersects $L_3$ if and only if $\beta = f(\alpha)$.
Since a rational fraction is fully determined by the image of $3$ points, you can find it by looking at $3$ particular lines intersecting $L_1,L_2,L_3$, for example the line intersecting $L_2$ and $L_3$ that is parallel to $L_1$ (that involves $\alpha = \infty$), the one intersecting $L_1$ and $L_3$ that is parallel to $L_2$ (that involves $\beta = \infty$), and the one that intersects $L_1$ and $L_2$ and is parallel to $L_3$
And if you have a fourth line $L_4$ then there will be another degree $1$ rational fraction $g$ such that $(IJ)$ intersects $L_4$ if and only if $\beta = g(\alpha)$.
So if there is a line that goes through $I$ and intersects $L_2,L_3,L_4$ then it has to go through $J = J_0+ \beta J_1$ with $\beta = f(\alpha) = g(\alpha)$, and this gives an equation on $\alpha$.
If $f(\alpha) = (a\alpha+b)/(c\alpha+d)$ and $g(\alpha) = (a'\alpha+b')/(c'\alpha+d')$, setting $f(\alpha) = g(\alpha)$ will give you the degree $2$ equation $(a\alpha+b)(c'\alpha+d') = (a'\alpha+b')(c\alpha+d)$.
A: To find the line, you need to set up the linear system of equations
$$\vec X_0 + \eta_1 \vec X_1 = \vec I_0 + \alpha \vec I_1$$
$$\vec X_0 + \eta_2 \vec X_1 = \vec J_0 + \beta \vec J_1$$
$$\vec X_0 + \eta_3 \vec X_1 = \vec K_0 + \gamma \vec K_1$$
$$\vec X_0 + \eta_4 \vec X_1 = \vec L_0 + \kappa \vec L_1$$
which imposes that the line $\vec X$ has a common point with all the other lines. This is a nonlinear system of $12$ equations (each equality has to hold for $3$ dimensions) in $13$ unknowns ($3$ values for $\vec X_0$, $3$ values for $\vec X_1$, and $7$ parameters $\eta_2,\ldots,\eta_4,\alpha,\beta,\gamma,\kappa$; there is a scaling factor in $\vec X_1$ that is solved by setting $\eta_1 = 1$). This underdetermined system may not have a solution, so i believe that not all sets of 4 lines have a line intersecting them all.
Edit:
You can further reduce the system by solving for $\vec X_0$ in the first equation,
$$ \vec X_0 = \vec I_0 + \alpha \vec I_1 - \vec X_1,$$
then replacing this in the remaining equations
$$\vec I_0 + \alpha \vec I_1 - \vec X_1 + \eta_2 \vec X_1 = \vec J_0 + \beta \vec J_1$$
$$ \vec I_0 + \alpha \vec I_1 - \vec X_1+ \eta_3 \vec X_1 = \vec K_0 + \gamma \vec K_1$$
$$\vec I_0 + \alpha \vec I_1 - \vec X_1+ \eta_4 \vec X_1 = \vec L_0 + \kappa \vec L_1$$
Further still:
Could this be useful? Six remaining unknowns with ? equations.
$$\frac{\vec J_0 + \beta \vec J_1}{\eta_2} = \frac{\vec K_0 + \gamma \vec K_1}{\eta_3} = \frac{\vec L_0 + \kappa \vec L_1}{\eta_4}$$
Or perhaps rearranging each and using each dimension:
$$\alpha \vec I_1  - \beta \vec J_1 + (1 - \eta_2) \vec X_1 = \vec J_0 - \vec I_0$$
$$
    \begin{bmatrix}
    {I_1}_x & {J_1}_x & {X_1}_x \\
    {I_1}_y & {J_1}_y & {X_1}_y \\
    {I_1}_z & {J_1}_z & {X_1}_z \\
    \end{bmatrix}
    \begin{Bmatrix}
    \alpha \\
    -\beta \\
    1 - \eta_2 \\
    \end{Bmatrix}
    =    
    \begin{Bmatrix}
    {J_0}_x - {I_0}_x \\
    {J_0}_y - {I_0}_y \\
    {J_0}_z - {I_0}_z \\
    \end{Bmatrix}
$$
...
However this is confusing as it suggests there is a single solution for each value which we know not to be the case. Alternately it suggests there are three solutions for each value instead of the expected two.
