The classical proof via differential geometry (which I might not have time to enter all of as I'm waiting for my flight; if not, I'll finish it when I get home) goes like this:
Suppose that the surface in $\mathbb{R}^3$ is smooth and parametrize it locally in the form $X(s,t)$ where the two rulings are defined by holding either $s$ or $t$ constant. Then This is a local argument, so, for simplicity, I'll assume that the domain of $X$ is a rectangle in the $st$-plane. Of course, one can reparametrize in $s$ and/or $t$ separately, and this will turn out to be useful at some point in the calculation.
The two tangent vector fields $X_s$ and $X_t$ are linearly independent and are the tangents to the two rulings. Since $X_{ss}$ is the acceleration of the $t$-ruling, it follows that $X_{ss} = f X_s$ for some function $f$. Similarly $X_{tt} = g X_t$ for some function $g$. Moreover
Note that, sinceif one reparametrized, using $\bar s$ and $\bar t$ instead of $s$ and $t$, then one would have
$$
X_{\bar s} = \frac{ds}{d\bar s}\ X_s\quad\text{and}\quad
X_{\bar s\bar s} = \left(f + \frac{d^2s}{d\bar s^2}\left(\frac{ds}{d\bar s}\right)^{-1}\right)\ X_{\bar s},
$$
with similar formulae for $X_{\bar t}$ and $X_{\bar t\bar t}$. This will be useful below.
Since the surface isdoes not lie in a plane, it follows that $X_{st}$ cannot be a linear combination of $X_s$ and $X_t$ (otherwise, the plane spanned by $X_s$ and $X_t$ would be fixed, and the surface would lie in plane). This means that $X_s$, $X_t$, $X_{st}$ is a basis of $\mathbb{R}^3$, and, as such, we havethere are equations of the form
$$
\begin{pmatrix} dX_s& dX_t & dX_{st}
\end{pmatrix}
= \begin{pmatrix} X_s& X_t & X_{st}
\end{pmatrix}
\begin{pmatrix} f\ ds& 0 & f_t\ ds\\\\
0 & g\ dt & g_s\ dt\\\\
dt & ds & f\ ds + g\ dt
\end{pmatrix}
$$
(The equations for $dX_{st}$ follow since $(X_{st})_s = (f X_s)_t = f_t\ X_s + f\ X_{st}$, etc.) By comparing partials, or by using the structure equations above (i.e., expanding out the consequences of $d(d(X_{st}))=0$, etc.), one sees that $d(f\ ds + g\ dt) = 0$, so that there must exist a function $h$ such that $f = h_s$$f = 2h_s$ and $g = h_t$$g = 2h_t$. (The coefficient of $2$ avoids some fractions later.)
The equation now becomes
$$
\begin{pmatrix} dX_s& dX_t & dX_{st}
\end{pmatrix}
= \begin{pmatrix} X_s& X_t & X_{st}
\end{pmatrix}
\begin{pmatrix} h_s\ ds& 0 & h_{st}\ ds\\\\
0 & h_t\ dt & h_{st}\ dt\\\\
dt & ds & h_s\ ds + h_t\ dt
\end{pmatrix}
$$$$
\begin{pmatrix} dX_s& dX_t & dX_{st}
\end{pmatrix}
= \begin{pmatrix} X_s& X_t & X_{st}
\end{pmatrix}
\begin{pmatrix} 2h_s\ ds& 0 & 2h_{st}\ ds\\\\
0 & 2h_t\ dt & 2h_{st}\ dt\\\\
dt & ds & 2h_s\ ds + 2h_t\ dt
\end{pmatrix}
$$
Moreover, the structure equations now imply that $d(e^hh_{st})=0$$d(e^{-2h}h_{st})=0$, so $h_{st} = ce^{-h}$$h_{st} = C\ e^{2h}$ for some constant $c$$C$. By adding a constant to $h$, one can reduce to the case that $C$ is one of $0$, $1$, or $-1$.
Consider the case $C=0$ (OOPS! Gotwhich needs to gobe treated separately in any case). Then $h_{st}=0$, so that, in particular, $f=2h_s$ is a function of $s$ alone and $g=2h_t$ is a function of $t$ alone. Using the change of variables formulae mentioned above, one can then change variables in $s$ so as to arrange that $f = 0$ and change variables in $t$ to arrange that $g=0$. Thus, the equations have reduced to
$$
d(X_s) = X_{st}\ dt,\qquad d(X_t) = X_{st}\ ds,\qquad d(X_{st})=0.
$$
Thus $X_{st} = v_3$ where $v_3$ is a constant vector. Then $d(X_s-tv_3) = 0$ and $d(X_t - s v_3) = 0$, so there exist constant vectors $v_1$ and $v_2$ so that $X_s = v_1 + t v_3$ and so that $X_t = v_2 + s v_3$. Finally, this implies that
$$
dX = X_s\ ds + X_t\ dt = (v_1+tv_3)\ ds + (v_2 + s v_3)\ dt
= d\bigl(sv_1+t v_2 + st v_3 \bigr),
$$
so that
$$
X = v_0 + sv_1+t v_2 + st v_3
$$
for some constant vector $v_0$. Thus, $X(s,t)$ parametrizes a hyperbolic paraboloid.
Now, consider the case $C\not=0$. The structure equations have become
$$
\begin{pmatrix} dX_s& dX_t & dX_{st}
\end{pmatrix}
= \begin{pmatrix} X_s& X_t & X_{st}
\end{pmatrix}
\begin{pmatrix} 2h_s\ ds& 0 & 2C\ e^{2h}\ ds\\\\
0 & 2h_t\ dt & 2C\ e^{2h}\ dt\\\\
dt & ds & 2\ dh
\end{pmatrix}
$$
where $h_{st} = C\ e^{2h}$. Moreover, one calculates
$$
d(e^{-2h} X_{st}) = 2C\ (X_s\ ds + X_t\ dt) = 2C\ dX
$$,
so it follows that $e^{-2h}X_{st} = 2C X + v_0$ for my flight! I'llsome constant vector $v_0$. In particular, showing that the vector-valued function $E_3 = e^{-2h}X_{st}$
takes values in a hyperboloid of $1$-sheet will finish the proof.
To this tonight when I get homeend, consider the new frame field
$$
\begin{pmatrix}E_1 & E_2 & E_3\end{pmatrix}
= \begin{pmatrix}e^{-h}X_{s} & e^{-h}X_{t} & e^{-2h}X_{st}\end{pmatrix}.
$$
Calculation shows that it satisfies the structure equation
$$
\begin{pmatrix} dE_1& dE_2 & dE_3
\end{pmatrix}
= \begin{pmatrix} E_1& E_2 & E_3
\end{pmatrix}
\begin{pmatrix} h_s\ ds-h_t\ dt& 0 & 2C\ e^h\ ds\\\\
0 & h_t\ dt - h_s\ ds & 2C\ e^h\ dt\\\\
e^h\ dt & e^h\ ds & 0
\end{pmatrix}.
$$
Note that the $3$-by-$3$ matrix of $1$-forms on the right takes values in the vector space
$\frak{g}$ consisting of matrices of the form
$$
\begin{pmatrix} x_1& 0 & 2C\ x_2\\\\
0 & -x_1 & 2C\ x_3\\\\
x_3 & x_2 & 0
\end{pmatrix}.
$$
This is, of course, the Lie algebra of the subgroup $O(Q)\subset GL(3,\mathbb{R})$ consisting of the matrices that satisfy $A^TQA = Q$ where
$$
Q = \begin{pmatrix} 0& 1 & 0\\\\
1 & 0 & 0\\\\
0 & 0 & -2C
\end{pmatrix}
$$
It follows that there is an invertible linear transformation $L$ of $\mathbb{R}^3$ such that the matrix $LE$ takes values in $O(Q)$, where $E = (E_1\ E_2\ E_3)$.) In particular $L$ carries the image of $E_3$ into the hyperboloid of $1$-sheet $2x_1x_2 - 2C x_3^2 = -2C$. By the remark above, it follows that $X(s,t)$ must be the image of this quadric under an affine transformation of $\mathbb{R}^3$, as was to be shown.
