最优化

本人总结的常用优化算法，主要讲解原理：http://download.csdn.net/detail/zk_j1994/9880349

主要包括梯度下降，共轭梯度法；牛顿法，拟牛顿法；信赖域方法，罚函数法。

# -*- coding: utf-8 -*-"""author: UniqueZ_file:   gradient descent algorithm, CG methoddate:   2017-06-23remark: 原方程为(x1)**2 + 2*(x2)**2"""import numpy as npimport matplotlib.pyplot as pltclass Gradient_method:    """ 使用的方向和梯度有关的优化算法: 梯度下降法 & 共轭梯度法 """    def __init__(self, x0, gradient_threshold, method):        """ 1. 初始化        x0: 初始迭代点        x: 当前迭代点        gradient_threshold: 算法终止的梯度阈值        method: "gradient" & "CG"        """        self.x = x0        self.gradient_threshold = gradient_threshold        self.method = method            def compute_original_fun(self, point):        """ 2. 计算点point的函数值        point: 一个点, 但不一定是当前的迭代点        """        value = point[0,0]**2 + 2*point[1,0]**2        return value            def compute_gradient(self):        """ 3. 计算当前迭代点的梯度值         value:当前迭代点的梯度值,  一个矩阵        """        value = np.mat([[0],[0]], np.float)        value[0,0] = 2*(self.x[0])        value[1,0] = 2*(self.x[1])        return value        def is_termination(self):        """ 4. 检验算法是否应该终止 """        if np.linalg.norm(self.compute_gradient()) < self.gradient_threshold:            return True        else:            return False        def compute_direction(self, d_last=None, g_last=None):        """ 5. 计算搜索方向         d_last: 上一次的迭代点的搜索方向        g_last: 上一次的迭代点的梯度        g: 当前迭代点的梯度        value: 当前迭代点的搜索方向, 一个矩阵        """        if self.method == "gradient":            value = - self.compute_gradient()        elif self.method == "CG":            g = self.compute_gradient()            value = -g + ((g.T*g)/(g_last.T*g_last))[0,0]*d_last        else:            raise(Exception("No method founded."))        return value        def compute_step(self, alpha = 1, row = 1/4):        """ 6. 计算步长         x: 当前迭代点, 一个行向量        alpha: 初始步长, 返回步长, 步长迭代过程中每次缩短1/2, 最多迭代20次, alpha变为初始alpha的百万分之一        row: 确定步长简单搜索法的一个参数        """        gradient = self.compute_gradient()        direction = - gradient        max_iter = 20        for i in range(max_iter):            if self.compute_original_fun(self.x + alpha*direction) <= self.compute_original_fun(self.x) + alpha*row*gradient.T*direction:                break            alpha = 1/3 * alpha            if alpha < 0.1:                alpha = 0.1            if i == max_iter - 1:                print("No valid step founded.")        return alpha        def compute_next_point(self, alpha, direction):        """ 7. 计算下一个迭代点 """        self.x = self.x + alpha * direction            def draw_result(self, point):        """ 8. 画出迭代点的变化过程 """        plt.figure()        plt.plot(range(len(point)), point, "y")        plt.title("min value's change")        return plt        def main(self, iter_point, max_iter=1000):        """ main        method: 梯度优化算法方面的方法: gradient(梯度下降) & CG(共轭梯度法)        max_iter: 最大迭代次数        x: 最优解        """        iter_point.append(np.linalg.norm(self.x))        for i in range(max_iter):            # 算法是否终止            if self.is_termination():                break                            # 计算方向direction            if self.method == "gradient":                direction = self.compute_direction()            elif self.method == "CG":                if i == 0:                    gradient = self.compute_gradient()                    direction = - gradient                else:                    direction = self.compute_direction(direction, gradient)                    gradient = self.compute_gradient()            else:                raise(Exception("No method founded."))                            # 计算步长alpha            alpha = self.compute_step()                        # 计算下一个迭代点x            self.compute_next_point(alpha, direction)            iter_point.append(np.linalg.norm(self.x))        return iter_point            if __name__ == "__main__":    # init    params = {"x0": np.mat([[100],[1.99]]), "gradient_threshold": 1e-6, "method": "gradient"}    gradient = Gradient_method(**params)        # comput     iter_point = []    iter_point = gradient.main(iter_point, max_iter=1000)        # draw    gradient.draw_result(iter_point).show()