最近在做Learning from data的习题，要解决一个线性优化问题，但是碰到一个坑，这里记录一下，具体题目可以参考Learning from data的第三章的Problem3.6，题目的大概意思是从线性规划角度解决分类问题，我先产生1000组数据。

import numpy as np
import matplotlib.pyplot as plt

#参数
rad=10
thk=5
sep=5

#n为产生点的个数,x1,y1为上半个圆环的坐标
def generatedata(rad,thk,sep,n,x1=0,y1=0):
    #上半个圆的圆心
    X1=x1
    Y1=y1

    #下半个圆的圆心
    X2=X1+rad+thk/2
    Y2=Y1-sep
    
    #上半个圆环的点
    top=[]
    #下半个圆环的点
    bottom=[]
    
    #后面要用到的参数
    r1=rad+thk
    r2=rad
    
    cnt=1
    while(cnt<=n):
        #产生均匀分布的点
        x=np.random.uniform(-r1,r1)
        y=np.random.uniform(-r1,r1)
        
        d=x**2+y**2
        if(d>=r2**2 and d<=r1**2):
            if (y>0):
                top.append([X1+x,Y1+y])
                cnt+=1
            else:
                bottom.append([X2+x,Y2+y])
                cnt+=1
        else:
            continue

    return top,bottom

n=1000
top,bottom=generatedata(rad,thk,sep,n)

X1=[i[0] for i in top]
Y1=[i[1] for i in top]

X2=[i[0] for i in bottom]
Y2=[i[1] for i in bottom]

plt.scatter(X1,Y1,s=1)
plt.scatter(X2,Y2,s=1)
plt.show()

显然这组数据是可分的，但是我带入scipy的linprog函数确始终解不出来。

from scipy.optimize import linprog

#加上偏移项
x1=[[1]+i for i in top]
x2=[[1]+i for i in bottom]

c=np.array([1.0,1.0,1.0])
A=np.concatenate((np.array(x1),np.array(x2)*(-1)))*(-1)
b=np.ones(n)*(-1.0)

result=linprog(c,A_ub=A,b_ub=b,options=dict(bland=True))
result
     fun: 71.923242295321842
 message: 'Optimization failed. Unable to find a feasible starting point.'
     nit: 1000
  status: 2
 success: False
       x: nan

看到有这样一个报错‘Optimization failed. Unable to find a feasible starting point.’，意思是没有解，但是显然有解啊，我就非常郁闷，上网查阅了一下，看到有人说100组数据以上linprog会解不出，于是我试了下50组数据。

n=50
top,bottom=generatedata(rad,thk,sep,n)
#加上偏移项
x1=[[1]+i for i in top]
x2=[[1]+i for i in bottom]

c=np.array([1.0,1.0,1.0])
A=np.concatenate((np.array(x1),np.array(x2)*(-1)))*(-1)
b=np.ones(n)*(-1.0)

result=linprog(c,A_ub=A,b_ub=b)
result
     fun: 0.68978181507748448
 message: 'Optimization terminated successfully.'
     nit: 55
   slack: array([  8.25582251,   1.39579263,   5.45099168,   7.89511373,
         4.07001065,   7.5398225 ,   6.36340605,   5.19580442,
         0.        ,   7.88915867,   7.21006349,   5.61077845,
         0.80448989,   5.84426297,   1.22033368,   5.00423168,
         4.2392862 ,   2.43914755,   0.67341515,   9.26841422,
         6.22546989,   8.54254805,   6.22844178,  12.44853646,
         8.24740109,   9.51136402,   4.8966531 ,  11.47498782,
         9.10895613,  11.81051557,  12.51911944,   7.09788684,
         3.43465829,  10.05898861,  11.83502706,  10.34493975,
         8.90589327,  11.95815476,   9.4623471 ,   2.78260225,
         6.77083725,   9.6965408 ,  10.10390994,   5.0151723 ,
         3.2516008 ,   4.95376907,  11.83342563,   7.88599283,
         2.53952582,   8.34889399])
  status: 0
 success: True
       x: array([ 0.        ,  0.        ,  0.68978182])

一看还真是这个问题，于是网上查阅了下有没有其他的工具，看到cvxopt是一个比较不错的凸优化工具，于是下载了后尝试了一下，效果不错，需要注意的一点是cvxopt的数据格式一定是浮点数，我一开始因为这个还报错了。

from cvxopt import matrix, solvers

n=1000
top,bottom=generatedata(rad,thk,sep,n)
#加上偏移项
x1=[[1]+i for i in top]
x2=[[1]+i for i in bottom]

c=np.array([1.0,1.0,1.0])
A=np.concatenate((np.array(x1),np.array(x2)*(-1)))*(-1)
b=np.ones(n)*(-1.0)

#要转换为matrix
c=matrix(c)
A=matrix(A)
b=matrix(b)

sol = solvers.lp(c,A,b)

print(sol)
     pcost       dcost       gap    pres   dres   k/t
 0:  3.4242e-01  1.0017e+03  2e+03  2e+00  7e+03  1e+00
 1:  3.4061e+01  4.7700e+04  8e+06  9e+01  3e+05  2e+02
 2:  7.0763e-01  4.3232e+02  6e+02  7e-01  3e+03  5e+01
 3:  8.4793e-01  5.2743e+02  1e+03  9e-01  3e+03  6e+01
 4:  1.8868e+00  7.3869e+01  1e+02  7e-02  3e+02  3e+01
 5:  1.8724e+00  3.3409e+00  5e+00  2e-03  8e+00  3e-01
 6:  1.4143e+00  1.8996e+00  1e+00  7e-04  3e+00  1e-01
 7:  1.4221e+00  1.8726e+00  1e+00  6e-04  2e+00  1e-01
 8:  1.3924e+00  1.5287e+00  3e-01  2e-04  8e-01  2e-02
 9:  1.3752e+00  1.3953e+00  5e-02  3e-05  1e-01  2e-03
10:  1.3731e+00  1.3733e+00  6e-04  4e-07  1e-03  3e-05
11:  1.3730e+00  1.3730e+00  6e-06  4e-09  1e-05  3e-07
12:  1.3730e+00  1.3730e+00  6e-08  4e-11  1e-07  3e-09
13:  1.3730e+00  1.3730e+00  6e-10  4e-13  1e-09  3e-11
Optimal solution found.
{'x': <3x1 matrix, tc='d'>, 'y': <0x1 matrix, tc='d'>, 's': <1000x1 matrix, tc='d'>, 'z': <1000x1 matrix, tc='d'>, 'status': 'optimal', 'gap': 6.385522247857282e-10, 'relative gap': 4.65067372803741e-10, 'primal objective': 1.3730316552293447, 'dual objective': 1.3730316554698367, 'primal infeasibility': 3.958299300636158e-13, 'dual infeasibility': 1.4801406159663262e-09, 'primal slack': 9.474745058411444e-13, 'dual slack': 9.915311930562724e-14, 'residual as primal infeasibility certificate': None, 'residual as dual infeasibility certificate': None, 'iterations': 13}

最后的结果是以字典形式存储的，比较重要的是’x’键所指的内容，我们的最优解对应的向量，这里作图看一下效果。

w=list((sol['x']))
#作图
r=2*(rad+thk)
X3=[-r,r]
Y3=[-(w[0]+w[1]*i)/w[2] for i in X3]

plt.scatter(X1,Y1,s=1)
plt.scatter(X2,Y2,s=1)
plt.plot(X3,Y3)
plt.show()

效果还是很不错的，以上就是这次报错的总结，所以以后大量数据的优化使用cvxopt比较好。