L1 optimization for fitting a curves

Question

I've done the Least square method of ellipse fitting. But I want a fitting method with less point. and convex optimization may be a strong one. From my study: $ax^2+bxy+cy^2+dx+ey+f=0$,

so $(a\ b \ c \ d \ e\ f)*(x^2 \ xy \ y^2 \ x \ y \ 1)=0, st.b^2-4*ac=1$, simplified as $\theta*z=0,st. \theta^TC\theta=1$, C = [0 0 2;0 -1 0;2 0 0], and the others is zeros.

Then, convert to unconstrained optimization $min_{\theta} \ |\theta z| + \lambda(\theta^TC\theta-1) $, I search goole, find a method called split Bregman. The equation becomes $min_{d,\theta} \ |d| + \lambda(\theta^TC\theta-1) +\frac{\mu}{2}|| d - \theta z -b ||_2^2, st. d=\theta z$

It's splited two part. One is $\theta^{k+1} = min_{\theta} \ \lambda (\theta^TC\theta-1)+\frac{\mu}{2}||d^k-\theta z-b^k||_2^2$, setting the derivative to zeros, the solution is $ \theta^{k+1} = \frac{\mu \theta^T(d^k-b^k)}{2\lambda C+\mu \theta^T \theta } $.

The other is $d^{k+1}=min_d \ |d^k| + \frac{\mu }{2} ||d^k - \theta^{k+1}z-b^k||_2^2$, according to the papaer, $d_j^{k+1} = shrink(\theta_j^{k}z + b_j^k,\frac{1}{\mu})$, and shrink operator is defined as $shrink(x,r)=\frac{x}{|x|}max(|x|-r,0)$.

The rest is loop. I write it as matlab, and find it produce hyperbolic, even I gives enough point.

I've checked the dimension error with exhaustion method.Also, I write the derivation so many times(I think it has no obvious errors). None of the code gives a ellipse. Below is my code with no running error.

Can you give me any advice? Thanks!

Refs: *Goldstein and Osher, The split Bregman method for L1 regularized problems SIAM Journal on Imaging Sciences 2(2) 2009

enter preformatted text here

clear all
clc

t = 0:0.1:4;
x = 3*sin(2*t+pi/4);
y = 2*cos(2*t);
f1 = [x;y]';
D = [f1(:,1).^2 f1(:,1).*f1(:,2) f1(:,2).^2 f1(:,1) f1(:,2) ones(size(f1,1),1)];

c0 = [0 0 2;0 -1 0;2 0 0];
c = [c0 zeros(3);zeros(3) zeros(3)];

n_max = 300;lambda = 1000;b=randn(size(D,1),1);
d=2*randn(size(D,1),1);
mu=1;
ii=1;
alpha_k0 = zeros(6,1);
alpha_k1 = pinv(lambda*c + mu*D'*D)*mu*D'*(d-b);

while abs(sum(D*alpha_k1)) > abs(sum(D*alpha_k0)) && (ii<n_max)
   alpha_k0 = alpha_k1;
   d = shrink(D,alpha_k1,b,mu);
   b = b - d + D*alpha_k1;
   alpha_k1 = pinv(lambda*c + mu*D'*D)*mu*D'*(d-b);
   ii = ii+1;
end

function out=shrink(D,alpha,b,mu)
   x = abs(D*alpha + b);
   out = max(abs(x)-1./mu,0).*(x./(abs(x)));
   out(isnan(out))=0;
end

enter preformatted text here