前言
在我上一篇博文《散点图下基于切比雪夫（Chebyshev）近似准则拟合直线》，介绍了如何用三点极小化最大残差的方法去拟合直线。该直线满足切比雪夫准则。本篇文章主要介绍如何用线性规划去拟合直线，该结果也满足切比雪夫准则，同时也证明了我上一篇博文的正确性。
线性规划
给定散点对(X,Y)，假设拟合的直线 ,
满足切比雪夫准则。则存在最大残差
，使得
。于是每一个点，我们有方程组：
给它取个负号，即两边同时乘以-1：
很明显，有多少个点，就有多少的方程组对。这是个线性规划问题，而我们的目的就是要求此线性规划使得成立的解。
松弛变量
此时，我们给方程组加上一个松弛变量，让其不等式变成等式：
其中就是松弛变量(relaxation variables)，且
。
现在，我们有变量。其中
是我们要求的。同时，它也说明，我们至少需要三个方程来求解。
代数法求解
如上图，假设有三个不等式方程构成如图的线性规划问题：
这里我们引入松弛变量：
当松弛变量取0时，我们发现满足方程的(x,y)在直线上。而由上述的不等式，我们知道在上图的灰色区域的点，必定是松弛变量满足
。而若两两松弛变量取0，其方程组求解出来的解就是其交点。而线性规划的目的，就是求解这些交点中的最优解。
现在，我们来求解方程组：
的最优解使得
求解前，我们先列出前提条件：。
情况一：
很明显，存在小于0的情况。
情况二：
若集合Y小于0，则。否则
，则
。很明显，
无法完全保证大于等于0。
情况三：
解方程组我们可以得到，b和有无穷多个解。同理
也一样。
情况四:
解得，
。
情况五:
解得，
，同时解得的其他松弛变量必须满足
情况六:
解得，
， 同时解得的其他松弛变量必须满足
情况七:
或者
解方程组我们可以得到，b和有无穷多个解。
情况八:
。
解得，
，
， 同时解得的其他松弛变量必须满足
情况九: 
解得，
，
， 同时解得的其他松弛变量必须满足
要求 ，我们只要比对以上每种情况中的最小者。
仔细看情况八和情况九，我们会发现，情况八和九就是我们上一篇所说的“三点极小化最大残差”。
因为表明直线平行于
的连线。而
将使得b的取值使得直线经过
构成的三角形的中线。
现我们选取情况五六七八进行比较，贴上代码：
function [ A,B,r ] = LinearProgram( X,Y )
[m,n]=size(X);
A=0;
B=0;
r=0;
A1=0;
B1=0;
r1=10000;
for i=1:m
for j=i+1:m
rp=(X(i,1)*Y(j,1)-X(j,1)*Y(i,1))/(X(i,1)-X(j,1));
ap=(Y(i,1)-((X(i,1)*Y(j,1)-X(j,1)*Y(i,1))/(X(i,1)-X(j,1))))/X(i,1);
if rp>=0 && ap>=0
if rp<r1
right=1;
for k=1:m
z0=Y(k,1)-ap*X(k,1)-rp;
z1=ap*X(k,1)-rp-Y(k,1);
if z0>0 || z1>0
right=0;
break;
end
end
if right==1
r1=rp;
A1=ap;
end
end
end
end
end
if r1~=0
A=A1;
B=B1;
r=r1;
end
A2=0;
B2=0;
r2=10000;
for i=1:m
for j=i+1:m
rp=(X(j,1)*Y(i,1)-X(i,1)*Y(j,1))/(X(i,1)+X(j,1));
ap=(Y(i,1)-((X(j,1)*Y(i,1)-X(i,1)*Y(j,1))/(X(i,1)+X(j,1))))/X(i,1);
if rp>=0 && ap>=0
if rp<r2
right=1;
for k=1:m
z0=Y(k,1)-ap*X(k,1)-rp;
z1=ap*X(k,1)-rp-Y(k,1);
if z0>0 || z1>0
right=0;
break;
end
end
if right==1
r2=rp;
A2=ap;
end
end
end
end
end
if r2~=0 && r>r2
A=A2;
B=B2;
r=r2;
end
A3=0;
B3=0;
r3=10000;
for i=1:m
for j=i+1:m
for k=1:m
if k~=j && k~=j
ap=(Y(i,1)-Y(j,1))/(X(i,1)-X(j,1));
rp=(Y(j,1)-Y(k,1)-ap*(X(j,1)-X(k,1)))*0.5;
b=Y(i,1)-ap*X(i,1)-rp;
if rp<r3
right=1;
for g=1:m
z0= Y(g,1)-ap*X(g,1)-rp-b;
z1=ap*X(g,1)-rp-Y(g,1)+b;
if z0>0 || z1>0
right=0;
break;
end
end
if right==1
A3=ap;
B3=b;
r3=rp;
end
end
end
end
end
end
if r3~=0 && r>r3
A=A3;
B=B3;
r=r3;
end
A4=0;
B4=0;
r4=10000;
for i=1:m
for j=i+1:m
for k=1:m
if k~=j && k~=j
ap=(Y(i,1)-Y(j,1))/(X(i,1)-X(j,1));
rp=(ap*(X(j,1)-X(k,1))-(Y(j,1)-Y(k,1)))*0.5;
b=rp-ap*X(i,1)+Y(i,1);
if rp<r4
right=1;
for g=1:m
z0= Y(g,1)-ap*X(g,1)-rp-b;
z1=ap*X(g,1)-rp-Y(g,1)+b;
if z0>0 || z1>0
right=0;
break;
end
end
if right==1
A4=ap;
B4=b;
r4=rp;
end
end
end
end
end
end
if r4~=0 && r>r4
A=A4;
B=B4;
r=r4;
end
end
第一张图为上篇文章拟合的结果，第二张为以上算法拟合的结果。