梯度下降法拟合球面

今天遇到了球面拟合的问题，正好在上coursera上的机器学习的课程，干脆就用梯度下降拟合一个球面试试吧

%% clear
clear; 
clc;
close all;

%% generate data 
theta = 0:2*pi/100:2*pi;
beta = pi/6:2*pi/100:5*pi/6;
r = 1;

x = r*cos(beta)'*cos(theta) + 1;
y = r*cos(beta)'*sin(theta) + 2;
z = r*sin(beta)'*ones(1,length(theta)) + 3;

%% add noise
n = randn(length(beta),length(theta));
n = n/max(max(n)) * 0.05;

z = z+n;

figure
mesh(x,y,z)

%% construct cost function
% J(a,b,c,R) = 1/(2*length(theta)*length(beta)) * Sigma((x-a)^2+(y-b)^2+(z-c)^2-R^2)^2
a = x(1);
b = y(1);
c = z(1);
R = 1.1;
alpha = 0.35;
N = 100;
J = zeros(1,N);

x = reshape(x, length(beta)*length(theta), 1);
y = reshape(y, length(beta)*length(theta), 1);
z = reshape(z, length(beta)*length(theta), 1);

tic
% iterate
for i=1:N
    deltaA = alpha * 1/(length(beta)*length(theta)) * sum(((x-a).^2+(y-b).^2+(z-c).^2-R^2).*(x-a)*-2);
    deltaB = alpha * 1/(length(beta)*length(theta)) * sum(((x-a).^2+(y-b).^2+(z-c).^2-R^2).*(y-b)*-2);
    deltaC = alpha * 1/(length(beta)*length(theta)) * sum(((x-a).^2+(y-b).^2+(z-c).^2-R^2).*(z-c)*-2);
    deltaR = alpha * 1/(length(beta)*length(theta)) * sum(((x-a).^2+(y-b).^2+(z-c).^2-R^2)*R*-2);
    
    a = a-deltaA;
    b = b-deltaB;
    c = c-deltaC;
    R = R - deltaR;
    
    J(i) = 1/(2*length(beta)*length(theta)) * sum(((x-a).^2 + (y-b).^2 + (z-c).^2-R^2).^2);
end

toc

figure
plot(J)

效果也还可以，100次迭代耗时0.044s，圆心x，y轴坐标误差很小，轴坐标差0.04, 半径差0.04。