3.10 拉格朗日乘子法 张满满 .mdl

MachineLearningMathBasic/3.10 拉格朗日算法基础 张满满.md

@@ -0,0 +1,218 @@
+#机器学习教学基础
+
+## 3.10 拉格朗日算法基础
+
+
+
+#3.10拉格朗日乘子法  张满满
+
+##讲义和规范
+
+•https://www.matongxue.com/madocs/939.html
+
+
+
+
+
+这篇博文中直观上讲解了拉格朗日乘子法和 KKT 条件，对偶问题等内容。 首先从无约束的优化问题讲起，一般就是要使一个表达式取到最小值： 
+
+
+
+
+
+*m**i**n**f*(*x*)minf(x)
+
+
+ 如果问题是 *m**a**x**f*(*x*)maxf(x)  ，这个是一个习惯。对于这类问题在高中就学过怎么做。只要对它的每一个变量求导，然后让偏导为零，解方程组就行了。 
+
+![二维线性可分示例图](https://images2018.cnblogs.com/blog/1189614/201804/1189614-20180409222801202-1630842805.png)
+
+     所以在极值点处一定满足 *d**f*(*x*)*d**x*=0df(x)dx=0  （只是必要条件，比如 *f*(*x*)=*x*3f(x)=x3  在 *x*=0x=0   处就不是极值点），然后对它进行求解，再代入验证是否真的是极值点就行了。对于有些问题可以直接通过这种方法求出解析解（如最小二乘法）。
+
+     但是也有很多问题解不出来或者很难解，所以就需要梯度下降法、牛顿法、坐标下降法之类的数值迭代算法了（感知机 、logistic 回归中用到）。
+
+     对于这些迭代算法就像下面这张图一样，我们希望找到其中的最小值。一个比较直观的想法是先找一个起点，然后不断向最低点靠近。就先把一个小球放到一个碗里一样。
+
+![迭代算法](https://images2018.cnblogs.com/blog/1189614/201804/1189614-20180409222844005-1090929262.bmp)
+
+     一开始要找一个起始点，然后确定走的方向和距离，最后还要知道什么时候停止。这三步中最难的应该是确定走的方向。走的慢点还可以接受，要是方向错了就找不到最小值了~。所以走的距离可以简单的设为一个比较小的值。起始点可以随机选一个 (*x*0,*y*0)(x0,y0)  。关键是方向，可以选择 (*x*0,*y*0)(x0,y0)   处的梯度的反方向，这是函数在这个点下降最快的方向（原因可以看知乎中忆臻的回答）。它是一个向量，然后它的大小就是走的距离，为了防止太大而走过头，导致不断在最小值附近震荡，需要乘上一个比较小的值（称为学习率），最终的停止条件就是梯度的大小很接近于 0（在极值点处的梯度大小就是 0）就行了。这种方法依靠梯度确定下降方向的方法叫做梯度下降法。
+
+ 对 *f*(*x*)f(x)   求极小值的流程就是：
+
+1. 随机选定*x*0x0  
+2. 得到函数在 *x*0x0   的梯度，然后从 *x*0x0    向前走一步。计算式是：*x**n**e**w*0=*x**o**l**d*0−*α*∇*f*(*x*)x0new=x0old−α∇f(x)  
+3. 重复第 2 步，直到梯度接近于 *o*o  （小于一个事先设定的很小的数），或者到达指定的迭代上限。
+
+![梯度下降法](https://images2018.cnblogs.com/blog/1189614/201804/1189614-20180409222921051-1966980537.bmp)
+
+     除了这种方法之外，其中第 2 步还可以这样做，固定 *x*x  , 把它作为常数。就变成只有一个变量的优化问题了，直接求导为 0 就可以得到最优点，向前走到 (*x*0,*y*1)(x0,y1)  处，然后固定 *y*1y1  , 对 *x*x   进行相同的操作。这种每次只优化一个变量的方法叫做坐标下降法。
+
+![坐标下降法](https://images2018.cnblogs.com/blog/1189614/201804/1189614-20180409222942020-1292610452.bmp)
+
+     然后就是进一步的，我们可能要在满足一定约束条件的情况下最小化目标函数，比如有一个等式约束： 
+
+
+
+
+
+*m**i**n**f*(*x*)*s*.*t*.*h*(*x*)=0minf(x)s.t.h(x)=0
+
+
+     解决这个的时候问题不能先用上面的方法求出 *f*(*x*)f(x)   的。因为这个问题的解可能根本不是 *m**i**n**f*(*x*)minf(x)
+
+
+
+![带约束的极值](https://images2018.cnblogs.com/blog/1189614/201804/1189614-20180409223000822-1093207370.png)
+
+     如图，其中的圆圈是指目标函数 *f*(*x*，*y*)f(x，y)   投影在平面上的等值线，表示在同一个圆圈上，黑线是约束条件*h*(*x*)=0h(x)=0  的函数图像。所以等值线与函数图像相交的点其实就是所有满足约束的点。那么极值点只有可能在等值线与函数图像相切的地方取到，因为如果在相交的地方取到，那么沿着 *h*(*x*)h(x)   的图像往前走或者往后走，一定还有其它的等值线与它相交，也就是 *f*(*x*,*y*)f(x,y)   的值还能变大和变小，所以交点不是极值点，只有相切的时候才有可能是极值点(不可能同时变大和变小了)。在相切的地方 *h*(*x*)h(x)   的梯度和 *f*(*x*,*y*)f(x,y)   的梯度应该是在同一条直线上的。（这一点可以这么想，在切点处两个函数的梯度都与切平面垂直，所以在一条直线上） 
+
+ 所以满足条件的极值点一定满足：∇*f*(*x*,*y*)=*λ*∇*h*(*x*,*y*)(*λ*=0是允许的，表示*f*(*x*,*y*)本身的极值点刚好在切点上)∇f(x,y)=λ∇h(x,y)(λ=0是允许的，表示f(x,y)本身的极值点刚好在切点上)  然后和原来的等式方程 *h*(*x*,*y*)=0h(x,y)=0   联立，然后只要解出这个方程组，就可以得到问题的解析解了。当然也存在解不出来的情况，就需要用罚函数法之类的方法求数值解了。
+
+     为了方便和好记，就把原来的优化问题写成 *f*(*x*,*y*)+*λ**h*(*x*,*y*)f(x,y)+λh(x,y)  的形式，然后分别对 *λ*,*x*,*y*λ,x,y  求偏导，并且令偏导为 0 就行了，和之前得到的方程组是一样的。这种方法叫拉格朗日乘数法。    
+
+ 如果有多个等式约束怎么办呢，如下图：
+
+![多个约束的极值](https://images2018.cnblogs.com/blog/1189614/201804/1189614-20180409223020563-581635895.png)
+
+     这里的平面和球面分别代表了两个约束 *h*1(*x*)h1(x)   和 *h*2(*x*)h2(x)  ，那么这个问题的可行域就是它们相交的那个圆。这里蓝色箭头表示平面的梯度，黑色箭头表示球面的梯度，那么相交的圆的梯度就是它们的线性组合（只是直观上的），所以在极值点的地方目标函数的梯度和约束的梯度的线性组合在一条直线上。所以就满足：
+
+
+
+
+
+
+
+∇*f*(*x*)=*λ**μ**i*∇*h**i*(*x*)=∑*i*=12*λ**i*∇*h**i*(*x*)*h*1(*x*)=0*h*2(*x*)=0∇f(x)=λμi∇hi(x)=∑i=12λi∇hi(x)h1(x)=0h2(x)=0
+
+
+ 大于2个约束的情况也一样。为了好记，将原来的约束的问题写成 *L*(*x*,*λ*)=*f*(*x*)+∑*i*−1*n**λ**i*∇*h**i*(*x*)L(x,λ)=f(x)+∑i−1nλi∇hi(x)   求偏导，然后让它们为 0。结果像上面一样直接列方程组是一样的。这个可以看做是一种简记，或者是对偶问题，这个函数叫做拉格朗日函数。
+
+  再进一步，如果问题中既有等式约束，又有不等式约束怎么办呢？对于： 
+
+
+
+
+
+*m**i**n**f*(*x*)*s*.*t*.*h*(*x*)=0*g*(*x*)≤0minf(x)s.t.h(x)=0g(x)≤0
+
+
+ 当然也同样约定不等式是 ⩽⩽   只要取反就行了。对于这个问题先不看等式约束，对于不等式约束和目标函数的图：
+
+![不等式约束](https://images2018.cnblogs.com/blog/1189614/201804/1189614-20180409223040318-1100903285.png)
+
+ 阴影部分就是可行域，也就是说可行域从原来的一条线变成了一块区域。那么能取到极值点的地方可能有两种情况：
+
+- 还是在*h*(*x*)h(x)   和 等值线相切的地方
+- *f*(*x*)f(x)   的极值点本身就在可行域里面。
+
+ 因为如果不是相切，那么同样的，对任意一个在可行域中的点，如果在它附近往里走或者往外走，*f*(*x*)f(x)  一般都会变大或者变小，所以绝大部分点都不会是极值点。除非这个点刚好在交界处，且和等值线相切；或者这个点在可行域内部，但是本身就是 *f*(*x*)f(x)  的极值点。如下图（维基百科上的图~）：
+
+![不等式约束下的极值](https://images2018.cnblogs.com/blog/1189614/201804/1189614-20180409223309966-1127532199.png)
+
+ 对于第一种情况，不等式约束就变成等式约束了，对*f*(*x*)+*λ**h*(*x*)+*μ**g*(*x*)f(x)+λh(x)+μg(x)  
+
+ 用拉格朗日乘子法：
+
+
+
+
+
+
+
+∇(*x*)+*λ*∇*h*(*x*)+*μ*∇*g*(*x*)=0*h*(*x*)=0*g*(*x*)=0*μ*≥0∇(x)+λ∇h(x)+μ∇g(x)=0h(x)=0g(x)=0μ≥0
+
+
+ 这里需要解释一下，为什么不是 *μ*≠0μ≠0  。后面的约束比前面的更强。看“不等式约束”那个图，我们已经知道了问题中的可行域是在 *g*(*x*)≤0g(x)≤0  的梯度是指向大于 0 的一侧，也就是不是可行域的一侧。而求的问题是极小值，所以 *f*(*x*)f(x)  的梯度方向，所以对它没有约束，那么为什么 *μ*μ   上。
+
+ 对于第二种情况，不等式约束就相当于没有，对 *f*(*x*)+*λ**h*(*x*)f(x)+λh(x)   用拉格朗日乘子法： 
+$$
+\begin{equation} 
+\begin{aligned} 
+\nabla(x) +\lambda\nabla h(x)+= 0\ 
+h(x) = 0 \ 
+g(x) \le 0\
+
+\end{aligned} 
+\end{equation} 
+$$
+ 最好把两种情况用同一组方程表示出来。对比一下两个问题，不同的是第一种情况中有 $\mu \ge 0$且 
+
+*g*(*x*)=0g(x)=0  , 第二种情况 *μ*=0μ=0   且 *g*(*x*)≤0g(x)≤0   综合两种情况，可以写成 *μ**g*(*x*)=0μg(x)=0  且 *μ*≥0μ≥0   且 *g*(*x*)≤0g(x)≤0  
+
+
+
+
+
+
+
+∇(*x*)+*λ*∇*h*(*x*)+*μ*∇*g*(*x*)=0*μ**g*(*x*)=0*μ*≥0*h*(*x*)=0*g*(*x*)≤0∇(x)+λ∇h(x)+μ∇g(x)=0μg(x)=0μ≥0h(x)=0g(x)≤0
+
+
+ 这个就是 KKT 条件。它的含义是这个优化问题的极值点一定满足这组方程组。（不是极值点也可能会满足，但是不会存在某个极值点不满足的情况）它也是原来的优化问题取得极值的必要条件，解出来了极值点之后还是要代入验证的。但是因为约束比较多，情况比较复杂，KKT 条件并不是对于任何情况都是满足的。要满足 KKT 条件需要有一些规范性条件（Regularity conditions），就是要求约束条件的质量不能太差，常见的比如：
+
+1. LCQ：如果 *h*(*x*)h(x)  和 *g*(*x*)g(x)  都是形如 *A**x*+*b*Ax+b  的仿射函数，那么极值一定满足 KKT 条件。
+2. LICQ：起作用的 *g*(*x*)g(x)   函数（即*g*(*x*)g(x)  相当于等式约束的情况）和 *h*(*x*)h(x)   函数在极值点处的梯度要线性无关，那么极值一定满足 KKT 条件。
+3. Slater 条件：如果优化问题是个凸优化问题，且至少存在一个点满足 *h*(*x*)=0h(x)=0  和*g*(*x*)=0g(x)=0  ，极值一定满足 KKT 条件。并且满足强对偶性质。
+
+
+
+定义：
+
+假设有自变量x和y，给定约束条件g(x,y)=c，要求f(x,y)在约束g下的极值。
+
+我们可以画出f的等高线图，如下图。此时，约束g=c由于只有一个自由度，因此也是图中的一条曲线（红色曲线所示）。显然地，当约束曲线g=c与某一条等高线f=d1相切时，函数f取得极值。
+ 两曲线相切等价于两曲线在切点处拥有共线的法向量。因此可得函数f(x,y)与g(x,y)在切点处的梯度（gradient）成正比。
+ 于是我们便可以列出方程组求解切点的坐标(x,y)，进而得到函数f的极值
+
+![1566981185187](C:\Users\这不是你的电脑\AppData\Roaming\Typora\typora-user-images\1566981185187.png)
+
+
+
+
+
+下面的例子是一道约束非线性优化问题，摘自高等数学规划作业题，原题为： 
+![这里写图片描述](http://www.57lai.com/images/IqqIqqIqqIqq/53735938/20161220124454514.jpg)
+
+
+
+```
+from scipy.optimize import fsolve
+import sympy
+espl=0.0001
+a0=-1
+def check(a,b,w):
+    s = a**2-b-3
+    print "a**2-b-3:",s
+    if s<float(w)/float(2):
+        return True
+    else:
+        return False
+def modify(w,v,x1,x2):
+    w = max(0,w-2*((x1-2)**2-x2))
+    v = v-2*(2*x1-x2-1)
+    return (w,v)
+def iterf(i,w,v,a):
+    b=(-w-v+2*a**2-6+4*a)/6
+    a=afs(w,v,b)
+    x1=a+2
+    x2=b+3
+    print "第%s次迭代：a=%s,b=%s,x1=%s,x2=%s" %(i,a,b,x1,x2)
+    wv = modify(w,v,x1,x2)
+    w = wv[0]
+    v = wv[1]
+    print "w=%s,v=%s" %(w,v)
+    #a = asolve(a,w,v,b)
+    return (w,v,a)
+def asolve(a,w,v,b):
+    #af = afs(a,w,v,b)
+    return fsolve(afs(a,w,v,b),a0)
+def afs(w,v,b):
+    res =  fsolve(lambda a: 2*a**3+(1-w-2*b-2)*a-v-2*b, 0)
+    #res = sympy.solve('a**3+((1-w)/2-b-1)*a-v/2-b')
+    flag = check(res,b,w)
+    print "是否满足：",flag
+    return res
+```
+
+
+