A continuación, utilizaré un método diferente al sugerido en su documento de referencia.
Como el problema de ejemplo indica que el elipsoide a identificar está centrado en el origen, la ecuación algebraica del elipsoide es
$ \mathbf{r}^T Q \mathbf{r} = 1 $
donde $\mathbf{r} = [x, y, z]^T $ y $Q$ es una simétrica positiva definida $3 \times 3$ matriz.
Dado un conjunto de puntos de datos $\{ P_i = (x_i, y_i, z_i) , i = 1, N\}, N \ge 6 $ se desea ajustar estos puntos al modelo
$ A x^2 + B y^2 + C z^2 + D xy + E xz + F yz + G = 0 $
Para ello, se construye una función de error $\mathbf{E}$ definido del siguiente modo
$ \mathbf{E} = \displaystyle \sum_{i=1}^N (A x^2 + B y^2 + C z^2 + D xy + E xz + F yz + G)^2 $
Y minimiza esta función. Para asegurarte de que obtienes un elipsoide válido con los parámetros $A$ a través de $G$ no todos cero, se hace la minimización condicionada a la restricción $A + B + C = 1 $
Sustituir por $A$ en la expresión anterior se obtiene
$ \mathbf{E} = \displaystyle \sum_{i=1}^N ( x_i^2 + B (y_i^2 - x_i^2) + C (z_i^2 - x_i^2) + D x_i y_i + E x_i z_i + F y_i z_i + G )^2 $
Diferenciación $\mathbf{E}$ con respecto a $B, C, D, E, F, G$ resulta en las siguientes ecuaciones para el mínimo
$ \displaystyle \sum_{i=1}^N ( x_i^2 + B (y_i^2 - x_i^2) + C (z_i^2 - x_i^2) + D x_i y_i + E x_i z_i + F y_i z_i + G ) (y_i^2 - x_i^2) = 0 $
$ \displaystyle \sum_{i=1}^N ( x_i^2 + B (y_i^2 - x_i^2) + C (z_i^2 - x_i^2) + D x_i y_i + E x_i z_i + F y_i z_i + G ) (z_i^2 - x_i^2) = 0 $
$ \displaystyle \sum_{i=1}^N ( x_i^2 + B (y_i^2 - x_i^2) + C (z_i^2 - x_i^2) + D x_i y_i + E x_i z_i + F y_i z_i + G ) (x_i y_i) = 0 $
$ \displaystyle \sum_{i=1}^N ( x_i^2 + B (y_i^2 - x_i^2) + C (z_i^2 - x_i^2) + D x_i y_i + E x_i z_i + F y_i z_i + G ) (x_i z_i) = 0 $
$ \displaystyle \sum_{i=1}^N ( x_i^2 + B (y_i^2 - x_i^2) + C (z_i^2 - x_i^2) + D x_i y_i + E x_i z_i + F y_i z_i + G ) (y_i z_i) = 0 $
$ \displaystyle \sum_{i=1}^N ( ( x_i^2 + B (y_i^2 - x_i^2) + C (z_i^2 - x_i^2) + D x_i y_i + E x_i z_i + F y_i z_i + G ) = 0 $
Lo anterior $6$ conducen al siguiente sistema lineal para encontrar $B, C, D, E, F, G$
Definir el vector $X = [B, C, D, E, F, G]^T $ las ecuaciones normales son
$ \mathbf{A} X = Y $
donde
$ \mathbf{A} = \displaystyle \sum_{i=1}^N A_i $
con
$A_i = \begin{bmatrix} (y_i^2 - x_i^2)^2 && (z_i^2 - x_i^2)(y_i^2 - x_i^2) && x_i y_i (y_i^2 - x_i^2) && x_i z_i (y_i^2 - x_i^2) && y_i z_i (y_i^2 - x_i^2) && (y_i^2 - x_i^2 ) \\ (y_i^2 - x_i^2) (z_i^2 - x_i^2) && (z_i^2 - x_i^2)^2 && x_i y_i (z_i^2 - x_i^2) && x_i z_i (z_i^2 - x_i^2) && y_i z_i (z_i^2 - x_i^2) && (z_i^2 - x_i^2 )\\ (y_i^2 - x_i^2) (x_i y_i) && (z_i^2 - x_i^2) (x_i y_i) && (x_i y_i)^2 && x_i^2 y_i z_i && x_i y_i^2 z_i && x_i y_i \\ (y_i^2 - x_i^2) (x_i z_i) && (z_i^2 - x_i^2) (x_i z_i) && x_i^2 y_i z_i && x_i^2 z_i^2 && x_i y_i z_i^2 && x_i z_i \\ (y_i^2 - x_i^2) (y_i z_i) && (z_i^2 - x_i^2) (y_i z_i) && (x_i y_i)(y_i z_i) && (x_i z_i) (y_i z_i) && (y_i z_i)^2 && y_i z_i \\ (y_i^2 - x_i^2) && (z_i^2 - x_i^2) && (x_i y_i) && (x_i z_i) && (y_i z_i) && 1 \end{bmatrix}$
Y
$ Y = \begin{bmatrix} - \sum x_i^2 (y_i^2 - x_i^2) \\ - \sum x_i^2 (z_i^2 - x_i^2) \\ - \sum x_i^2 (x_i y_i) \\ - \sum x_i^2 (x_i z_i) \\ - \sum x_i^2 (y_i z_i) \\ - \sum x_i^2 \end{bmatrix}$
Ahora tenemos nuestro sistema lineal que se puede resolver para el vector parámetro $X$ utilizando la eliminación de Gauss-Jordan, que es una rutina estándar. Habiendo encontrado $X$ se conocen todos los parámetros del elipsoide y se puede escribir en forma cuadrática como
$ r^T Q r = 0 $
donde $r = [x, y, z, 1]^T$ y
$ Q = \begin{bmatrix} A && D/2 && E/2 && 0 \\ D/2 && B && F/2 && 0 \\ E/2 && F/2 && C && 0 \\ 0&&0&&0 && G \end{bmatrix} $
Ahora defina el $3 \times 3$ matriz
$ Q_0 = \begin{bmatrix} A && D / 2 && E / 2 \\ D / 2 && B && F/2 \\ E/2 && F/2 && C \end{bmatrix}$
Lo que nos permite escribir la ecuación del elipsoide ajustado en la forma
$ \mathbf{r}^T Q_0 \mathbf{r} = - G $
donde $p = [x, y, z]^T $
Dividiendo ambos lados por $(-G ) $ nos da
$ \mathbf{r}^T Q_1 \mathbf{r} = 1$
donde $Q_1 = Q_0 / ( - G) $
Utilizando los datos de la primera tabla que has proporcionado, he obtenido
$Q_1 = \begin{bmatrix} 0.25 && -9 \times 10^{-7} && -2.7 \times 10^{-6} \\ -9\times 10^{-7} && 0.999999 && 1.56 \times 10^{-5} \\ -2.7 \times 10^{-6} && 1.56 \times 10^{-5} && 4 \end{bmatrix} $