Least squares fitting of an ellipse

least squaresMATLAB

I'm trying to use least squares to fit an ellipse to some data. I can follow the steps outlined in:

http://www.scielo.br/scielo.php?script=sci_arttext&pid=S1982-21702015000200329

Up to equation 12, however I don't understand how the author can solve the least squares problem in Matlab as per equation 13 from the paper, below.

$$ v=\begin{bmatrix}
x^{2} & y^{2} & z^{2} & 2xy & 2xz & 2yz & 2x & 2y & 2z
\end{bmatrix}\backslash\text{ones(n)} $$

If as per the previous document we write the equation to be solved as:

$$ \phi v = L $$

Where L is length n containing 1's, I assume as it should be a unit ellipse with magnitude 1.

Rearranging to solve gives:

$$ v = (\Phi\Phi^{T})^{-1}\Phi^{T}L $$

The Matlab mldivide (backslash) operator is equivalent to writing:

$$ A^{-1}b=A\backslash b $$

However I can't see how the equation is solved written in this manner in Matlab, would you not need to solve as per:

D = [ x .* x, ...
      y .* y, ...
      z .* z, ...
  2 * x .* y, ...
  2 * x .* z, ...
  2 * y .* z, ...
  2 * x, ...
  2 * y, ... 
  2 * z ];

% Solve the normal system of equations.
v = ( D' * D ) \ ( D' * ones( size( x, 1 ), 1 ) );

Would someone be able to point out the gap in my understanding?

Best Answer

According to elements of SLT P.64, given a linear relationship between $\mathbf{X}$ and $\mathbf{y}$ represented as $\mathbf{y} = \mathbf{X}^{T}\beta + \epsilon$ has a closed form unique solution $\hat{\beta} = (\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{y}$ iff $\mathbf{X}$ is invertible. It seems that this formulation is quite close to the one in eq. 14 $\hat{v} = (\mathbf{\Omega}^{T}\mathbf{P}\mathbf{\Omega})^{-1}\mathbf{\Omega}^{T}\mathbf{P}\ell$. The only difference is that $\mathbf{P}$ is a weight matrix giving emphasis on particular solutions over other existing ones. In eq. 22 it shows that is solving a linear equation of the form $\mathbf{J}\mathbf{d} = \mathbf{h}$ using projected coordinate descent. If my understanding is correct this is an approximate solution and not an exact one. To answer your question more conretely about eq. 13, it seems like a vector of dim $1\times 9$ solved with an identity matrix $\mathbf{I}_{9\times 9}$, now I don't have matlab but maybe matlab turns the vector into a matrix using broadcasting? But, after eq. 10 the authors describes a $n\times u$ design matrix which is missing in eq. 13. That could be a mistake although I haven't gone through great detail over the paper? Maybe not exactly the solution you were looking for but maybe some of the hints can help you move along?

Related Question