拉格朗日松弛方法

For a better view of math equations,
go to https://www.zybuluo.com/kailaix/note/490410.

拉格朗日松弛方法的基本思想是将“顽固”的约束（通常是等式约束）整合到目标函数中去。如果原来目标函数是一个线性函数，等式约束是一个线性约束，拉格朗日方法要求我们能够轻松求解给定区域下的线性目标函数问题。

拉格朗日方法的理论基础是凸优化中的对偶原理，需要注意的是，即使最后取得 max min 时候的极值点(argmin argmax $L(\mu)$)不是原来函数的极值点，但是在一定条件下，我们能得到一样的函数值。

另一个需要注意的是，对于去掉线性约束后约束区域无界的情况，需要小心处理（见以下笔记）；虽然在数学上容易处理，但是在实践中很容易造成溢出。

下面的代码是上面两个问题的求解，注意这里我们是手动求最小值，也可以借助现有优化包例如scipy.optimize。

from scipy import optimize as opt
import numpy as np
from math import pi,cos,sin
import cPickle


class ILRSolver(object):
    def __init__(self,f,marker):
        self.f = f
        self.stepsize = 1.0
        self.cnt = 0
        self.iter_cnt = 0
        self.mu = 0.0
        self.marker = marker
        self.L = []

    def _fmin(self):
        if self.marker == 1:
            if self.mu >= 1:
                self.x = np.array([0,0,0]);return 
            else:
                self.x = np.array([1,1,1]);return 
        if self.marker == 2:
            if self.mu < 1:
                self.x = np.array([1,1,1]);return 
            elif self.mu<2:
                self.x = np.array([0,1,1]);return 
            elif self.mu<3:
                self.x = np.array([0,0,1]);return 
            else:
                self.x = np.array([0,0,0]);return 

    def update_coef(self):
        delta = self.x.sum() - 1.5
        if delta > 0:
            self.mu += self.stepsize
        elif delta < 0:
            self.mu -= self.stepsize


    def update(self):
        self.update_coef()
        self._fmin()
        self.fx_old = self.fx
        self.fx = self.f(self.x)

        if self.fx == self.fx_old:
            self.cnt += 1

        if self.cnt >= 10 or self.iter_cnt % 10==0:
            self.cnt = 0
            self.stepsize /= 2.0

    def _nonstop(self):
        print "---------------------\n Current infomation:"
        print "iter: "+ str(self.iter_cnt)
        print "stepsize: "+str(self.stepsize)
        print "x: "+ str(self.x)
        print "mu: "+str(self.mu)
        print "f(x): "+ str(self.fx) 
        print "g(x): " + str(self.x.sum() - 1.5)
        obj_func = lambda x:self.f(x)+self.mu*(x[0]+x[1]+x[2]-1.5)
        print "objf(x): "+str(obj_func(self.x))
        self.L.append(obj_func(self.x))
        self.iter_cnt += 1
        # raw_input('press to continue')
        if self.stepsize <= 1e-5:
            print 'STEPSIZE TOO SMALL'
            return False
        if self.x[0]+self.x[1]+self.x[2] == 1.5:
            print 'OPTIMAL'
            return False
        return True

    def iterate(self):
        self._fmin()
        self.fx = self.f(self.x)
        while self._nonstop():
            self.update()

    def save_data(self):
        with open(str(self.marker)+'.dat','wb') as fp:
            cPickle.dump(self.L, fp)

f1 = lambda x: - ( x[0]**2 + x[1]**2 + x[2]**2 )
f2 = lambda x: - (x[0]+2*x[1]+3*x[2])

ilr1 = ILRSolver(f1,1)
ilr2 = ILRSolver(f2,2)

ilr1.iterate()
ilr1.save_data()

ilr2.iterate()
ilr2.save_data()

下面是两个对偶函数最大值的收敛形态。