线性规划下的直线拟合

2021-09-10 22:04:16 阅读：368 来源： 互联网

前言

在我上一篇博文《散点图下基于切比雪夫（Chebyshev）近似准则拟合直线》，介绍了如何用三点极小化最大残差的方法去拟合直线。该直线满足切比雪夫准则。本篇文章主要介绍如何用线性规划去拟合直线，该结果也满足切比雪夫准则，同时也证明了我上一篇博文的正确性。

线性规划

给定散点对(X,Y)，假设拟合的直线 $y_{i}=ax_{i}+b$ , $i=0,1,2,3,...$ 满足切比雪夫准则。则存在最大残差 $r_{max}>0$ ，使得 $r_{max}-\left | y_{i}-ax_{i}-b \right |\geqslant 0$ 。于是每一个点，我们有方程组：

$\left\{\begin{matrix} r_{max}-(y_{i}-ax_{i}-b)\geqslant 0\\ r_{max}+(y_{i}-ax_{i}-b)\geqslant 0 \end{matrix}\right.$

给它取个负号，即两边同时乘以-1：

$\left\{\begin{matrix} y_{i}-ax_{i}-b-r_{max}\leqslant 0\\ ax_{i}-r_{max}-y_{i}+b\leqslant 0 \end{matrix}\right.$

很明显，有多少个点，就有多少的方程组对。这是个线性规划问题，而我们的目的就是要求此线性规划使得 $Min(r_{max})$ 成立的解。

松弛变量

此时，我们给方程组加上一个松弛变量，让其不等式变成等式：

$\left\{\begin{matrix} y_{i}-ax_{i}-b-r_{max}+z_{i0}=0\\ ax_{i}-r_{max}-y_{i}+b+z_{i1}=0 \end{matrix}\right.$

其中 $z_{i0},z_{i1}$ 就是松弛变量(relaxation variables)，且 $z_{i0},z_{i1}\geqslant 0$ 。

现在，我们有变量 $(a,b,r_{max}),(z_{00},z_{01},z_{10},z_{11},...z_{n0},z_{n1})$ 。其中 $(a,b,r_{max})$ 是我们要求的。同时，它也说明，我们至少需要三个方程来求解。

代数法求解

如上图，假设有三个不等式方程构成如图的线性规划问题：

$\left\{\begin{matrix} a_{1}x+b_{1}-y\leqslant 0\\ a_{2}x+b_{2}-y\leqslant 0\\ a_{3}x+b_{3}-y\leqslant 0 \end{matrix}\right.$

这里我们引入松弛变量：

$\left\{\begin{matrix} a_{1}x+b_{1}-y+z_{1}=0\\ a_{2}x+b_{2}-y+z_{2}=0\\ a_{3}x+b_{3}-y+z_{3}=0 \end{matrix}\right.$

当松弛变量取0时，我们发现满足方程的(x,y)在 $a_{i}x+b_{i}+y=0$ 直线上。而由上述的不等式，我们知道在上图的灰色区域的点，必定是松弛变量满足 $z_{i}\geqslant 0$ 。而若两两松弛变量取0，其方程组求解出来的解就是其交点。而线性规划的目的，就是求解这些交点中的最优解。

现在，我们来求解方程组： $\left\{\begin{matrix} y_{i}-ax_{i}-b-r_{max}+z_{i0}=0\\ ax_{i}-r_{max}-y_{i}+b+z_{i1}=0 \end{matrix}\right.$ $i=0,1,2,3,...$

的最优解使得 $Min(r_{max})$

求解前，我们先列出前提条件： $r_{max}\geqslant 0,z_{i}\geqslant 0,i=1,2,3,4,....$ 。

情况一： $a=0,b=0,r_{max}=0$

很明显， $z_{i}$ 存在小于0的情况。

情况二： $a=0,b=0,z_{i0}=0$

若集合Y小于0，则 $r_{max}<0$ 。否则 $r_{max}=Max(y_{i})$ ，则 $z_{i0}=Max(y)-y_{i},z_{i1}=Max(y)+y_{i}$ 。很明显， $z_{i1}$ 无法完全保证大于等于0。

情况三： $a=0,z_{i0}=0,z_{j0}=0$

解方程组我们可以得到，b和 $r_{max}$ 有无穷多个解。同理 $a=0,z_{i1}=0,z_{j1}=0$ 也一样。

情况四: $a=0,z_{i0}=0,z_{j0}=0$

解得 $r_{max}=Max(\frac{abs(y_{i}-y_{j})}{2})$ ， $b=\frac{(y_{i}+y_{j})}{2}$ 。

情况五: $b=0,z_{i0}=0,z_{j0}=0$

解得 $r_{max}=Max(\frac{x_{i}y_{j}-x_{j}y_{i}}{x_{i}-x_{j}})$ ， $a=\frac{(y_{i}-\frac{x_{i}y_{j}-x_{j}y_{i}}{x_{i}-x_{j}})}{x_{i}}$ ，同时解得的其他松弛变量必须满足 $z_{i}\geqslant 0$

情况六: $b=0,z_{i0}=0,z_{j1}=0$

解得 $r_{max}=Max(\frac{y_{i}x_{j}-y_{j}x_{i}}{x_{i}+y_{j}})$ ， $a=\frac{y_{i}-\frac{y_{i}x_{j}-y_{j}x_{i}}{x_{i}+x_{j}}}{x_{i}}$ ，同时解得的其他松弛变量必须满足 $z_{i}\geqslant 0$

情况七: $z_{i0}=0,z_{j0}=0,z_{k0}=0$ 或者 $z_{i1}=0,z_{j1}=0,z_{k1}=0$

解方程组我们可以得到，b和 $r_{max}$ 有无穷多个解。

情况八: $z_{i0}=0,z_{j0}=0,z_{k1}=0,i\neq j\neq k$ 。

解得 $a=\frac{y_{i}-y_{j}}{x_{i}-x_{j}}$ ， $r_{max}=\frac{y_{j}-y_{k}-a(x_{j}-x_{k})}{2}$ ， $b=y_{i}-ax_{i}-r_{max}$ ，同时解得的其他松弛变量必须满足 $z_{i}\geqslant 0$

情况九: $z_{i1}=0,z_{j1}=0,z_{k0}=0,i\neq j\neq k$

解得 $a=\frac{y_{i}-y_{j}}{x_{i}-x_{j}}$ ， $r_{max}=\frac{a(x_{j}-x_{k})-(y_{j}-y_{k})}{2}$ ， $b=r_{max}-ax_{i}+y_{i}$ ，同时解得的其他松弛变量必须满足 $z_{i}\geqslant 0$

要求 $Min(r_{max})$ ，我们只要比对以上每种情况中的最小者。

仔细看情况八和情况九，我们会发现，情况八和九就是我们上一篇所说的“三点极小化最大残差”。

因为 $a=\frac{y_{i}-y_{j}}{x_{i}-x_{j}}$ 表明直线平行于 $(x_{i},y_{i}),(x_{j},y_{j})$ 的连线。而 $Min(r_{max})$ 将使得b的取值使得直线经过 $(x_{i},y_{i}),(x_{j},y_{j}),(x_{k},y_{k})$ 构成的三角形的中线。

现我们选取情况五六七八进行比较，贴上代码：

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

第一张图为上篇文章拟合的结果，第二张为以上算法拟合的结果。

标签：直线,right,rp,线性规划,ap,松弛,拟合,end,解得
来源： https://blog.csdn.net/u010050735/article/details/120222752

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

ICode9

线性规划下的直线拟合

前言

线性规划

松弛变量

代数法求解

情况一：

情况二：

情况三：

情况四:

情况五:

情况六:

情况七: 或者

情况八:。

情况九: