I have used the method described below many times without issue to determine the major and minor axes of the "data ellipse". However, for this set of data (Fig. A) the method described below is failing to determine physical results. My question is, what is happening, why is the method giving the result it is?

enter image description here

Fig. B above is a plot of the data, in Cartesian coordinates. It is a plot of the square-root of the data $R(\alpha)=\sqrt k$ as a function of angle $\alpha$. The x-coordinate is $R\cos \alpha$ and the y-coordinate is $R\sin \alpha$. As seen, the data approximates the shape of an ellipse. Material heterogenieties from which this data is derived causes the scatter/deviation seen.

Method for determination of major and minor axes of the "data ellipse" using the principal of least squares: The directional data as a function of azimuthal angle $\alpha$ and principal values $k'_i$ (the eigenvalues of the data tensor) is given by Eqn(1) according to theory:

$$\tag{1} \frac{1}{k}=\frac{\cos^2 \alpha}{k'_1}+\frac{\sin^2 \alpha}{k'_2}$$

Eqn(1) is the equation of an ellipse if $R=\sqrt k$ is plotted vs the angle of measurement. Let $K_i$ equal $1/k_i$, where $k_i$ is the data measurements dependent on the angle of measurement $\alpha_i$. Then, $$\tag{2} K_i = K_{11}\cos ^2 \alpha_i + 2K_{12} \cos \alpha_i \sin \alpha_i + K_{22} \sin^2 \alpha_i$$

If $\rho_i$ equals the actual measured values of $1/k_i$ as a function of the angle $\alpha_i$, and if Eqn(2) is the ellipse fitted to these data, then according ot the principal of least squares,

$$\tag{3} G(K_{11}, K_{22}, K_{12})=\sum_i (\rho_i-K_i)^2$$

Taking the derivatives of Eqn(3) with respect to $K_{11}$, $K_{22}$, and $K_{12}$, and then some algebraic manipulation gives:

$$\tag{4} \begin{matrix} K_{11}\sum_i \cos^4 (\alpha_i)+K_{22}\sum_i \sin^2 (\alpha_i) \cos^2 (\alpha_i) =\sum_i \rho_i \cos^2 (\alpha_i) \\ K_{11}\sum_i \cos^2 (\alpha_i) \sin^2 (\alpha_i)+K_{22}\sum_i \sin^4 (\alpha_i) =\sum_i \rho_i \sin^2 (\alpha_i)\\ K_{12}\sum_i 2\cos^2 (\alpha_i) \sin^2 (\alpha_i) =\sum_i \rho_i \cos (\alpha_i) \sin (\alpha_i) \end{matrix}$$

To arrive at Eqn(4) we assumed data measured at equal intervals in 180 degrees, and therefore the sums of the terms with odd powers are zero.

Since the data given above were measured every 22.5 degrees, starting with 0 degree (ie $i = 1, 2, ..., 16$), the solution of Eqn (4) for $K_{11}, K_{22}, K_{12}$ is:

$$\tag{5} \begin{matrix} K_{11}=(3/16)\sum_i \rho_i \cos^2 (\alpha_i)-(1/16)\sum_i \rho_i \sin^2(\alpha_i)\\ K_{22}=(3/16)\sum_i \rho_i \sin^2 (\alpha_i)-(1/16)\sum_i \rho_i \cos^2(\alpha_i)\\ K_{12}=(1/4)\sum_i \rho_i \cos (\alpha_i) \sin (\alpha_i) \end{matrix}$$

With the values of $K_{11}, K_{22}, K_{12}$ now solved for using Eqn (5), we can solve for the normal and cross coefficients of the data's 2D tensor using:

$$\tag{6} [k]=[K^{-1}]=\frac{1}{det(K)}\cdot adj(K)$$


$$\tag{7} det(K)=K_{11}\cdot K_{22} - K^2_{12}$$


$$\tag{8} adj(K)=\begin{bmatrix}k_{22}& -k_{12}\\ -k_{21}&k_{11} \end{bmatrix}$$

The principal values of the 2D data tensor are found by solving the matrix eigenvalue problem; the eigenvalues $\lambda_i \equiv k'_i$ are given by

$$\tag{9} \lambda_i=1/2[k_{11}+k_{22}\pm [(k_{11}-k_{22})^2+4k^2_{12}]^{1/2}]$$

Below are the values I calculate for Eqns (5) - (9):

enter image description here

As seen, a principal value (one of the axes of the ellipse) of -10 md is being calculated which is physically not possible.

Now, if we estimate the angle of the direction of the maximum data value (the orientation of the major axis of the data ellipse), and then use non-linear regression to minimize the square of the errors in the measured values of $k$ directly, I find the data is better approximated with $k'_1 = 32$ and $k'_2=1.1$ (with the major axis orientation angle being $\theta_1 = 23$ degrees). Here is the 'best fit' ellipse drawn with these values on top of the data given above (for the data 22.5 through 337.5, every 45 degree):

enter image description here

As seen, this data is well approximated by the ellipse using the non-linear method. But I am puzzled as to why the least squares method is not working in this instance. Why is the least squares method giving the (erroneous) result it is? Is there a limitation of the least squares method as function of degree of anisotropy/heterogeneity of the data?

Edit 1 based on whuber's comment:

The following was taken from a textbook: $z_1, z_2,...,z_q$ are $q$ unknown quantities which, when multiplied by the coefficients $a_1, a_2, ..., a_q$, give the quantity $M$, according to the equation $$\tag{a} M=a_1z_1+a_2z_2+...+a_qz_q$$ With assumed values of the coefficients $a_1, a_2, ..., a_q$ we can set up an experiment and measure $M$. Owing to experimental errors we shall not measure the true value of $M$ but shall find a value $M+v$, where $v$ is the error. We could set up other experiments, $p$ in all, with different values of the $a$’s and write $$\tag{b} \begin{matrix} v_1=a_{11}z_1+a_{12}z_2+...+a_{1q}z_q-M_1 \\ v_2=a_{21}z_1+a_{22}z_2+...+a_{2q}z_q-M_2\\ ... \\ v_p=a_{p1}z_1+a_{p2}z_2+...+a_{pq}z_q-M_p \end{matrix}$$ If now $p>q$, there are a number of ways in which $q$ equations can be chosen from the $p$ equations and for each set we could, by assuming the errors to be zero, solve for the $q$ unknowns $z_1,z_2,…,z_q$, but the sets would not agree exactly. The principle of least squares states that the values of the $z$’s should be such that $v_1^2+v_2^2+⋯+v_p^2$ is a minimum. Equations (b) may be written $$\tag{c} v_i=a_{ij}z_j-M_i (i=1,...,p; j=1,...,q)$$ Or $$\tag{d} \bf v=az-M$$ where $\bf v, z$ and $\bf M$ are single-column matrices, and $\bf a$ is the matrix of the $a_{ij}$ coefficients. Since $v_i v_i$ is a minimum we have $$\tag{e} v_i \frac{\partial v_i}{\partial z_j}=0$$ Or, in view of Equation (c), $$\tag{f} v_i a_{ij}=0$$ Eliminating the $v_i$ between (c) and (f), and changing the dummy suffix from $j$ to $k$, we obtain $$\tag{g} (a_{ik}z_k-M_i)a_{ij}=0$$ or $$\tag{h} (a_t)_{ji}(a_{ik}z_k-M_i)=0$$ Which may be written without suffixes as $$\tag{i} \bf a_taz-a_tM=0$$ $\bf (a_ta)$ is a square matrix and has a reciprocal. Hence, multiplying before each term by $\bf (a_t a)^{-1}$, we solve for $\bf z$ and find $$\tag{j} \bf z=(a_t a)^{-1}a_tM$$

So, in the context of this thread’s question, the $z$'s are the unknown quantities we want to solve for. The coefficients $a_1,a_2,…,a_q$ are the components (or elements) of the vector-matrix-vector product of the direction cosines of the orientation of the measurement times the $k_{ij}$ matrix times the direction cosines matrix, e.g.,

$$\tag{k} M_1=(\cos \theta_1, \sin \theta_1, 0)\cdot \begin{pmatrix} k_{11} \ k_{12} \ 0 \\ k_{21} \ k_{22} \ 0 \\ 0 \ 0 \ k_{33} \end{pmatrix} \cdot \begin{pmatrix} \cos \theta_1 \\ \sin \theta_1 \\ 0 \end{pmatrix}$$

$$\tag{l} M_1=\cos^2 \theta_1 k_{11} + 2\cos \theta_1 \sin \theta_1 k_{12}+\sin^2 \theta_1 k_{22}$$

Therefore, the matrix $\bf a$, the matrix of the $a_{ij}$ coefficients, when using the trig identity $2 \cos \theta \sin \theta=2 \sin \theta$, is:

$$\tag{m} \begin{pmatrix} cos^2 \theta_1 \ 2 \sin \theta_1 \ sin^2 \theta_1 \\ ... \ ... \ ... \\ cos^2 \theta_n \ 2 \sin \theta_n \ sin^2 \theta_n \end{pmatrix}$$

For this case, using the angles given above for each measurement, $\bf a$ = enter image description here

Therefore, $\bf a_t a$ = enter image description here

and $\bf (a_t a)^{-1}$ = enter image description here

So we now have for $\bf (a_t a)^{-1}a_t$ = enter image description here

Finally, $\bf z$ is given by Eqn(j), we find: enter image description here

These are the same values that were calculated above using Equation (5).

  • $\begingroup$ Are you really measuring the $k_i = 1/R_i^2$?? Typically one measures the $x$ and $y$ coordinates of the points, or perhaps their distances and angles from an origin, or even multiple distances from reference points. None of those would lead to your formulation of the problem. This suggests the problem is a mismatch between your data and the model you are using. $\endgroup$
    – whuber
    Commented May 9, 2019 at 15:58
  • $\begingroup$ @whuber $R_i=\sqrt(k_i)$, so $k_i = R_i^2$ $\endgroup$
    – Armadillo
    Commented May 9, 2019 at 17:11
  • $\begingroup$ Could you explain how you measure $R_i^2$ directly? $\endgroup$
    – whuber
    Commented May 9, 2019 at 17:12
  • $\begingroup$ $k$ is an experimentally measured valued, $k = q\mu L/(A\Delta p)$. $q$ is volumetric flow rate, $\mu$ is fluid viscosity, $L$ is the length of flow, $A$ is cross-sectional area normal to flow, $\Delta p$ is the pressure drop driving flow. $k$ derives from a second rank tensor. Here we are looking at data within a plane (2D). The nature of the experiment used to calculate $k$ has the flow vector $\vec q$ in the direction of the measurement and the pressure gradient vector $\vec \Delta p$ at some oblique angle due to the anisotropy (cont.)... $\endgroup$
    – Armadillo
    Commented May 9, 2019 at 17:23
  • $\begingroup$ @whuber ...Therefore, what is really measured is $k = \vec q \mu L/(A |\vec \Delta p| \cos \theta)$ $\endgroup$
    – Armadillo
    Commented May 9, 2019 at 17:24


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.