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}
The key idea is that the matrix:
$$B=\begin{bmatrix}1&1&1\\1&1&0\\1&0&0\end{bmatrix}$$
transforms a vector in the base B to a vector in the standard basis and its inverse transform a vector in the standard basis to a vector in the base B:
$$v_S=B v_B \iff v_B=B^{-1}v_S$$
Thus, since the linear transformation in the standard basis is expressed by:
$$L=\begin{bmatrix}0&0&4\\3&5&-2\\1&1&4\end{bmatrix} \quad w_S=Lv_S$$
we have
$$w_S=Lv_S \implies Bw_B=LBv_B\implies w_B=B^{-1}LBv_B$$
and therefore the matrix
$$B^{-1}LB=\begin{bmatrix}1&1&1\\1&1&0\\1&0&0\end{bmatrix}^{-1}.\begin{bmatrix}0&0&4\\3&5&-2\\1&1&4\end{bmatrix}.\begin{bmatrix}1&1&1\\1&1&0\\1&0&0\end{bmatrix}$$
represents the linear transformation with respect to the new basis.
Best Answer
Fixed my issue - turns out I was misusing the formula.
It worked when I used $ Q' = (M^{-1})^T Q M^{-1} $ instead.