Be careful in how you try to solve this. Tools like fsolve will NOT be correct. In fact, fsolve will converge to the all zero solution in general, because that ALWAYS solves the problem exactly.
Even if you have 6 points for those 6 parameters (and there is no noise at all in your data) the general solution to the problem is A=B=C=D=E=F=0. See that this set of coefficients satisfies all 6 data points, yielding an exact result, regardless of the set of (x,y) pairs.
Now, suppose you have 60 points? Still the zero solution I have shown will always be exact. In fact, that will always be true that a zero solution exists.
Is there a valid way to solve the problem, yielding the conic form, when that is possible? Again, note my statement about when it is possible.
So lets try an example. Here, I'll create a simple circle in the plane, with center at the origin, and radius 2.
t = rand(6,1)*2*pi;
x = 2*cos(t);
y = 2*sin(t);
Now, create the matrix M. Assume column vectors for x and y.
M = [x.^2, x.*y, y.^2, x, y, ones(size(x))];
Does a non-zero solution exist for the linear homogeneous problem? Thus:
M*coefs = 0
A unique conic form will result only when the matrix M has rank exactly 5, so 1 less than the number of coefficients to be estimated. If the rank of M was 6, then no solution could ever exist beside the all zero solution. If the rank of M is less than 5, then there will be infinitely many conic forms possible, thus infinitely many sets of coefficients.
rank(M)
ans =
5
Great. The rank of M is perfect. A solution exists, and we can find it in theory. But how?
NS = null(M)
NS =
0.2357
3.0531e-16
0.2357
-1.138e-15
-2.1164e-15
-0.94281
That may not look that useful, but it gives us the solution. We need to recognize that the coefficients of a circle as you have written it are not in fact unique. We can multiply by any arbitrary constant and still have a valid circle. So, if we scale the resulting coefficients, to make the coefficient of x^2 be 1, then we see
NS/NS(1)
ans =
1
1.2953e-15
1
-4.828e-15
-8.979e-15
-4
(For a more fully problem, I might choose to scale it to make the coefficient with maximum absolute value 1. Scaling the problem here only serves to make it obvious that we have indeed solved for a circle, and gotten the one we started with.)
Thus
x^2 + y^2 -4 = 0
The other coefficients are down in the floating point trash.
Finally, suppose we change the problem slightly, by trashing one data point.
yhat = y;
yhat(1) = 2;
Mhat = [x.^2, x.*yhat, yhat.^2, x, yhat, ones(size(x))];
rank(Mhat)
ans =
6
As I said above, if the rank of M is 6, no solution is possible. Only the all zero solution exists.
This points out an issue with the above scheme. Suppose you have noise in your points? Any noise at all may be sufficient. Even a trivial amount of noise is a problem.
x = 2*cos(t) + randn(size(t))*1e-12;
y = 2*sin(t) + randn(size(t))*1e-12;
M = [x.^2, x.*y, y.^2, x, y, ones(size(x))];
rank(M)
ans =
6
So M is not in fact singular. But it is very close to that status.
cond(M)
ans =
7.8251e+13
You might arguably decide this is close enough to being numerically singular that you might blindly go ahead and use an SVD to solve the problem. Looking at the singular values of M, I would concede that M is indeed close to a singular matrix, of numerical rank very close to 5.
svd(M)
ans =
8.1681
5.0482
3.5868
2.039
0.20622
1.0438e-13
So...
[U,S,V] = svd(M);
NS = V(:,6);
NS/NS(1)
ans =
1
1.1607e-12
1
-4.2942e-13
4.7678e-12
-4
So we see success there.
By way of comparison, looking at the perturbed Mhat from before, we see:
svd(Mhat)
ans =
9.9209
5.0067
3.6948
0.9926
0.44018
0.07361
cond(Mhat)
ans =
134.78
Mhat clearly indicates no conic was ever to be gained from that data. The problems will arise when your noise is larger in magnitude.
t = rand(100,1)*2*pi;
x = 2*cos(t) + randn(size(t))*1e-2;
y = 2*sin(t) + randn(size(t))*1e-2;
M = [x.^2, x.*y, y.^2, x, y, ones(size(x))];
cond(M)
ans =
332.28
svd(M)
ans =
30.313
19.798
15.015
13.9
13.072
0.09123
So M is absolutely not singular, but I have lots of data here, though it is also moderately noisy. The requirement to even attempt a conic form solution is that the last singular value is MUCH smaller then the penultimate one. That is the case here, but not by much.
Even so, trying brute force, we will get:
[U,S,V] = svd(M);
NS = V(:,6);
NS/NS(1)
ans =
1
-0.00052009
0.99826
0.0042785
0.001349
-3.9975
Well, we almost recovered the original circle. The coefficients would indicate a nearly circular ellipse. We have B^2-4*A*C is less than zero, so the conic will be elliptical, although the major and minor axes are very close in magnitude.
Ok, is there any more to the above analysis? HELL YES! In fact, there is a lot more to say, but I will not be convinced to write the book it deserves.
If there is noise in your data points, all of this becomes now an error in variables problem. So if both x and y have noise as they were measured, it will introduce biases in the estimates of the parameters. This topic alone is worth a book. I'll leave it to say that noise in your data may cause problems, beyond those described above.
