First the equations for all vectors $x$ on line $g$ and all vectors $y$ on line $h$:
$$
\begin{align}
g: x &= a + \lambda b \quad \\
h: y &= c + \mu d
\end{align}
\quad (*)
$$
The difference vector between two of those vectors is
$$
D = y - x = c - a + \mu d - \lambda b
$$
The length of $D$ squared is:
\begin{align}
q(\lambda, \mu)
&= D \cdot D \\
&= (c - a)^2 + (\mu d-\lambda b)^2 + 2 (c-a)\cdot(\mu d - \lambda b) \\
&= (c - a)^2 + \mu^2 d^2 + \lambda^2 b^2 - 2 \mu \lambda (d\cdot b)
+ 2 \mu ((c-a)\cdot d) - 2 \lambda ((c-a)\cdot b) \\
\end{align}
The gradient of $q$ is:
$$
q_\lambda = 2\lambda b^2 - 2\mu (d\cdot b)-2((c-a)\cdot b) \\
q_\mu = 2\mu d^2 - 2\lambda (d\cdot b)+2((c-a)\cdot d)
$$
It should vanish for local extrema of $q$ which leads to the system
$$
A u = v
$$
with
$$
A =
\left(
\begin{matrix}
b^2 & - d\cdot b \\
- d \cdot b & d^2
\end{matrix}
\right)
\quad
u =
\left(
\begin{matrix}
\lambda \\
\mu
\end{matrix}
\right)
\quad
v =
\left(
\begin{matrix}
(c-a) \cdot b \\
-(c-a) \cdot d
\end{matrix}
\right)
$$
A unique solution exists for
$$
0 \ne \mbox{det A} = b^2 d^2 - (d \cdot b)^2
$$
Note that the dot operator stands for the scalar product.
That solution is
$$
u = A^{-1} v
$$
with
$$
A^{-1} = \frac{1}{\mbox{det } A}
\left(
\begin{matrix}
d^2 & d\cdot b \\
d \cdot b & b^2
\end{matrix}
\right)
$$
Inserting the found values for $\lambda$ and $\mu$ into equations $(*)$ will provide you the two vectors, whose points are closest to each other.
Example:
$$
a = (2,0,0), b=(1,1,1), c=(0,1,-1), d=(-1,0,-1)
$$
leads to
$$
\lambda = 1, \mu = -2.5, x_\min = (3,1,1), y_\min=(2.5,1,1.5), D=(-0.5,0,0.5)
$$
The image renders the $x$ values in green, the $y$ values in purple, and $D$ in red.