Finding the Transformation to a Canonical form for a Quadric Surface
It took a while, but I have an answer.
This answer describes the transformation of a hyperbolic paraboloid into canonical form. The procedure is similar for other quadric shapes.
An affine transform $\mathbf{T}$ is of the form: \begin{equation} \mathbf{T} = \begin{bmatrix} T_{11} & T_{12} & T_{13} & T_{14} \\ T_{21} & T_{22} & T_{23} & T_{24} \\ T_{31} & T_{32} & T_{33} & T_{34} \\ 0 & 0 & 0 & 1 \end{bmatrix} \end{equation}
Every quadric may be transformed via affine transformations into either a diagonal matrix or a near diagonal matrix, i.e. \begin{equation} \mathbf{X}^T \mathbf{T}^T \mathbf{A} \mathbf{T} \mathbf{X} = \mathbf{X}^T \mathbf{D} \mathbf{X} \end{equation} where $\mathbf{D}$ is a (near-)diagonal matrix. A near-diagonal matrix is symmetric and is of the form: \begin{equation} \mathbf{D} = \begin{bmatrix} D_{11} & 0 & 0 & 0 \\ 0 & D_{22} & 0 & 0 \\ 0 & 0 & 0 & D_{43} \\ 0 & 0 & D_{43} & 0 \end{bmatrix} \end{equation}
The canonical form of a hyperbolic paraboloid listed in Xu (Table 1, pp518) is not strictly in near-diagonal form: \begin{align} 0 &= z - xy \\ &= \mathbf{X}^T \mathbf{D'} \mathbf{X} \\ &= \mathbf{X}^T \begin{bmatrix} 0 & -0.5 & 0 & 0 \\ -0.5 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0.5 \\ 0 & 0 & 0.5 & 0 \end{bmatrix} \mathbf{X} \end{align} but may be re-written in near-diagonal form as: \begin{align} 0 &= -x^2 + y^2 + z \\ &= \mathbf{X}^T \mathbf{D} \mathbf{X} \\ &= \mathbf{X}^T \begin{bmatrix} -1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0.5 \\ 0 & 0 & 0.5 & 0 \end{bmatrix} \mathbf{X} \end{align}
The Xu form $\mathbf{D'}$ and the near-diagonal form $\mathbf{D}$ are related via an affine transform: \begin{align} \mathbf{D'} = \mathbf{M}^T \mathbf{D} \mathbf{M} \end{align} where: \begin{equation} \mathbf{M} = \begin{bmatrix} 0.5 & -0.5 & 0 & 0 \\ 0.5 & 0.5 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix} \end{equation}
Since the Xu's canonical form can be found from the near-diagonal form, our immediate task is to find the transformation $\mathbf{T}$ that casts a given hyperbolic paraboloid $\mathbf{A}$ into near-diagonal form. Although Dupont recommends determining the transform via Gaussian Reductions to avoid radicals, perhaps the easiest first step is to remove the rotation via an Eigendecomposition as suggested by Levin.
Consider the principle submatrix $\mathbf{A}_u$. Since it is symmetric, of rank 2 and has known signs for the eigenvalues, an eigendecomposition may be found such that \begin{align} \mathbf{A}_u = {\mathbf{R}_u}^T \mathbf{\Lambda} \mathbf{R}_u \end{align} where $\mathbf{R} \in SO(3)$ and $\mathbf{\Lambda}$ is in the form: \begin{equation} \mathbf{\Lambda} = \begin{bmatrix} \lambda_+ & 0 & 0 \\ 0 & \lambda_- & 0 \\ 0 & 0 & 0 \\ \end{bmatrix} \end{equation} where $\lambda_+ > 0$ and $\lambda_- < 0$. Therefore, the first affine transform may be built as: \begin{align} \mathbf{T}_1 = \begin{bmatrix} \mathbf{R}_u & \mathbf{0} \\ \mathbf{0} & 1 \end{bmatrix} \end{align}
The second transform forces the eigenvalues to unit magnitude and it can be verified that: \begin{align} \mathbf{T}_2 = \begin{bmatrix} \frac{1}{\sqrt{\lambda_+}} & 0 & 0 & 0 \\ 0 & \frac{1}{\sqrt{\|\lambda_-\|}} & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix} \end{align}
Let $\mathbf{A}_2 = {\mathbf{T}_2}^T {\mathbf{T}_1}^T \mathbf{A} \mathbf{T}_1 \mathbf{T}_2$, which will be of the form: \begin{equation} \mathbf{A}_2 = \begin{bmatrix} 1 & 0 & 0 & a \\ 0 & -1 & 0 & b \\ 0 & 0 & 0 & c \\ a & b & c & d \end{bmatrix} \end{equation}
It can be verified that the transform required to cast into near-diagonal form is: \begin{equation} \mathbf{T}_3 = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ \frac{-a}{c} & \frac{-a}{b} & \frac{-1}{2c} & \frac{-d}{2c} \\ 0 & 0 & 0 & 1 \end{bmatrix} \end{equation}
Therefore, the complete transformation required to cast the hyperbolic paraboloid into Xu's canonical form is: \begin{align} \mathbf{T} = \mathbf{T}_1 \mathbf{T}_2 \mathbf{T}_3 \mathbf{M} \end{align}