Our goal is to `complete the square' in 3D. Basically, this involves putting the quadratic form into the form of
\[
Q = \mathbf{x}^T A \mathbf{x},
\]

where $A$ is a real symmetric matrix. Thus, expand
\[
\begin{pmatrix}
x & y & z
\end{pmatrix}
\begin{pmatrix}
a_{11} & a_{12} & a_{13} \\
a_{12} & a_{22} & a_{23} \\
a_{13} & a_{23} & a_{33}
\end{pmatrix}
\begin{pmatrix}
x \\
y \\
z\end{pmatrix}
\]

and equate to the left hand side of the quadratic surface equation. This gives
\[
A = \begin{pmatrix}
5 & -5 & 1 \\
-5 & 11 & -5 \\
1 & -5 & 5
\end{pmatrix}.
\]

Now you can use the following result: the eigenvectors of $A$ give the directions of the principle axes of the surface (stationary values of its radius), and the eigenvalues give squares of the corresponding radii. In the case of $A$, we can verify that the eigenvalue and vector pairs are
\begin{gather*}
\lambda_1 = 16, \qquad v_1 = \frac{1}{\sqrt{6}} (1, -2, 1)^T \\
\lambda_2 = 4, \qquad v_2 = \frac{1}{\sqrt{2}} (-1, 0, 1)^T \\
\lambda_3 = 1, \qquad v_3 = \frac{1}{\sqrt{3}} (1, 1, 1)^T
\end{gather*}

We know from the theory of matrices that we can diagonalize the matrix $A$ in the form $A = S^{-1} D S$ where $S$ is a matrix of eigenvectors and $D$ is a matrix of eigenvalues. However, if $A$ is real and symmetric, then we can normalize all the eigenvectors, and this allows us to use $A = S^T D S$ due to the orthogonality of $S$. In our case,
\[
S = \begin{pmatrix}
\frac{1}{\sqrt{6}} & -\frac{2}{\sqrt{6}} & \frac{1}{\sqrt{6}} \\
\frac{-1}{\sqrt{2}} & 0 & \frac{1}{\sqrt{2}} \\
\frac{1}{\sqrt{3}} & \frac{1}{\sqrt{3}} & \frac{1}{\sqrt{3}}
\end{pmatrix}
\qquad
D = \begin{pmatrix}
16 & 0 & 0 \\
0 & 4 & 0 \\
0 & 0 & 1
\end{pmatrix}
\]

Now suppose that we'd like to perform an orthogonal transformation, $\mathbf{x} = S\mathbf{x}'$. We get for the quadratic form
\[
Q = \mathbf{x}^T A \mathbf{x} = \mathbf{x'}^T \underbrace{S^T A S}_{A'} \mathbf{x'}.
\]

However, $S^T A S = D$. So the quadratic form now takes a very simple form of
\[
Q = \mathbf{x'}^T D \mathbf{x'} = \lambda_1 x_1'^2 + \lambda_2 x_2'^2 + \lambda_3 x_3'^2
\]

so in the particular case of $Q = 1$,
\[
\left(\frac{x_1'}{1/2}\right)^2 + x_2'^2 + \left(\frac{x_3'^2}{2}\right)^2 = 1,
\]

so we conclude that the length of the semi-major lengths are $2, 1, 0.5$, and for the length corresponding to $2$, the eigenvector is along $(1, 1, 1)^T$.