Ayuda con el ajuste específico del elipsoide por mínimos cuadrados

Question

Siguiendo el método del artículo de Li y Griffiths al que se puede acceder libremente aquí He creado el siguiente problema de juguete para ver si estoy realizando los cálculos correctamente. A continuación describo el problema y doy mis cálculos intermedios y finales. Espero entender dónde y qué estoy haciendo mal.

El problema : Utilizando un elipsoide centrado en el origen definido por la ecuación $$\tag{1} 0.25x^2+1y^2+4z^2=1,$$ los siete puntos siguientes $p_i(x_i,y_i,z_i)$ que satisfacen la Ecuación (1):

A continuación se muestra un gráfico de estos puntos junto con el elipsoide:

Mi intento de ajustar un elipsoide a este conjunto de puntos siguiendo el método de Li y Griffiths : Con $k$ igual a 4, la matriz de restricciones $C_1$ es:

La matriz $DD^T$ (véase la sección 3 del documento, ecuación (11)):

La matriz $A=S_{11}-S_{12}S_{22}^{-1}S_{12}^T$ encontrada en la Ecuación(15) del documento, es:

La matriz $C_1^{-1}A$ es:

El mayor valor propio $\lambda$ (el único valor propio positivo) de la matriz $C_1^{-1}A$ y el vector propio asociado $u_1$ es:

Vector $u_2=-S_{22}^{-1}S_{12}^Tu_1$ es igual a 0,62. Así, con el vector de parámetros del elipsoide $u=(u_1,u_2)$ el elipsoide mejor ajustado es el siguiente $$\tag{2} ax^2+by^2+cz^2+2fyz+2gxz+2hxy+d=0$$ con los coeficientes:

El valor de $kJ-I^2$ donde $I$ y $J$ se definen en las ecuaciones (2) y (3) del documento, es 0,3 mientras que $k=4$ . El valor de $k$ necesarios para hacer $kJ-I^2=1$ es de 5,298, con el nuevo vector de parámetros asociado $u$ es

Esto es claramente erróneo. ¿En qué me he equivocado? ¿Qué he hecho mal?

Answer 1

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} $