用SVR拟合圆

一、 实验原理

SVR与SVM不同，SVR的样本点只有一类，它的目标是使得所有样本点到最优超平面的总偏差最小，而SVM的目标是让所有样本点到最优超平面的距离最大。
求圆的圆心和半径时，先把非线性问题转成线性问题：

因此本实验的优化目标为：

约束条件为：

其中ξ是松弛因子，ε为点到圆的距离偏差。
转换成拉格朗日函数为：

通过对拉格朗日核心因子（w，m，ξi+，ξi-）进行梯度下降，以及对拉格朗日乘子（αi+，αi-，βi+，βi-）进行梯度上升，即可得到最优的圆。

二、 实验过程

1. 读取数据

import os
path="数据集"
files=os.listdir(path)
coord_x=[]
coord_y=[]
for file in files:
    if not os.path.isdir(file):
        tmp_x=[]
        tmp_y=[]
        f = open(path+"/"+file)
        lines=f.readlines()
        for line in lines:
            coordinates=line.strip().split(" ")
            tmp_x.append(float(coordinates[0]))
            tmp_y.append(float(coordinates[1]))
        coord_x.append(tmp_x)
        coord_y.append(tmp_y)    

2. 数据可视化

import matplotlib.pyplot as plt
import numpy as np
def paint(x):
    fig,ax=plt.subplots(1,3)
    for i in range(0,3):
        plt.subplot(1,3,i+1)
        plt.scatter(coord_x[i+x],coord_y[i+x],s=1)
        plt.axis('square')
        plt.title("datasets"+str(i+x))
    fig.tight_layout()
    plt.show() 
paint(0)
paint(3)
paint(6)       

3. 用对偶下降法求解圆心和半径：

from tqdm import tqdm
import numpy as np
import math
import random

def dual_desc_ed2(total_epoch,X,Y,step,epsilon):
    length=len(X)

    # 随机生成拉格朗日因子
    alpha_plus=np.random.rand(length)
    beta_plus=np.random.rand(length)
    alpha_minus=np.random.rand(length)
    beta_minus=np.random.rand(length)

    z=[[0]*length for _ in range(3)]
    # 构造Z矩阵
    for i in range(length):
        z[0][i]=X[i]
        z[1][i]=Y[i]
        z[2][i]=z[0][i]**2+z[1][i]**2 
    z=np.array(z)    

    a=random.random()
    b=random.random()
    r=random.random()

    # 构造拉格朗日核心因子
    # reshape(-1，1)是把行向量转成列向量
    w=np.squeeze(np.random.random(3).reshape(-1,1))
    m=a**2+b**2-r**2
    xi_plus=np.random.rand(length)
    xi_minus=np.random.rand(length)

    # 设置超参数C
    C=random.random()
    # 开始迭代
    for _i in tqdm(range(total_epoch)):
        #　梯度下降求核心参数
        w-=step*(w+(alpha_plus*z).sum(axis=1)-(alpha_minus*z).sum(axis=1)) # 1.
        w[2]=1
        m-=step*(alpha_plus.sum()-alpha_minus.sum()) # 2.
        xi_plus-=step*(C-beta_plus-alpha_plus) # 3.
        xi_minus-=step*(C-beta_minus-alpha_minus) # 4.
        # 梯度上升求拉格朗日参数    
        beta_plus+=step*(-xi_plus) # 5.
        beta_minus+=step*(-xi_minus) # 6.
        alpha_plus+=step*(np.dot(w.T,z)+m-epsilon-xi_plus) # 7.
        alpha_minus+=step*(-epsilon-xi_minus-np.dot(w.T,z)-m) # 8.
           
    return w,m

radius=[]
centre_x=[]
centre_y=[]
total_epoch=100000
step=0.01
epsilon=0.01
# 求解圆心和半径
for i in range(9):
    result=dual_desc_ed2(total_epoch,coord_x[i],coord_y[i],step,epsilon)
    x=result[0][0]*-0.5
    y=result[0][1]*-0.5
    try:
        r=math.sqrt(pow(x,2)+pow(y,2)-result[1])
    except BaseException:
        print("radius of dataset {} errors".format(i))    
    centre_x.append(x)
    centre_y.append(y)
    radius.append(r)
    print("i={},x={},y={},r={}".format(i,x,y,r))              

4.绘制图形

# 绘制图形
from matplotlib.patches import Circle
def paintCircle(x):
    fig,ax=plt.subplots(1,3)
    for i in range(0,3):
        # print("(x,y)=({},{}),r={}".format(centre_x[i+x],centre_y[i+x],radius[i+x]))
        circle=Circle(xy=(centre_x[i+x],centre_y[i+x]),radius=radius[i+x],color='r',fill=False)
        plt.subplot(1,3,i+1)
        plt.scatter(coord_x[i+x],coord_y[i+x],s=1)
        plt.gca().add_artist(circle)
        plt.axis('square')
        plt.title("datasets"+str(i+x))
    fig.tight_layout()
    plt.show() 
paintCircle(0)
paintCircle(3)
paintCircle(6)    

三、 实验结果

迭代次数为100000时的效果如下，可以看到，epoch=100000时，对圆和3/4圆的拟合效果都比较好，但是对1/4圆的拟合效果则比较一般。

四、总结分析

1.错误分析：本次实验我把下图中数组Z的组成看错了，以为每一列为（x2，y2，x2+y2），实际上应该是（x，y，x2+y2），导致后续梯度计算时无法正确计算。

2.参数设置分析：
本次实验需要计算的参数包括两大类：拉格朗日核心因子（w，m，ξi+，ξi-）以及拉格朗日乘子（αi+，αi-，βi+，βi-），一共有八百多个，初始化时均选择用np.random.random()进行随机初始化。

此外，本实验还有3个超参数需要设置：学习率step，以及下图中的C和ε，除了学习率预设为0.01，C和ε我都用的是np.random.random()进行随机初始化。