最优化 - 牛顿优化算法

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

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

# -*- coding: utf-8 -*-
"""
author: UniqueZ_
file:   牛顿法, 基于DFP的拟牛顿法
date:   2017-06-24
remark: 原方程为(x1)**2 + 2*(x2)**2
"""
import numpy as np
import matplotlib.pyplot as plt

class Newton_method:
    """ 使用的方向和牛顿方向有关的优化算法: 牛顿法 & 基于DFP的拟牛顿法 """
    def __init__(self, x0, gradient_threshold, method, H=None):
        """ 1. 初始化
        x0: 初始迭代点
        x: 当前迭代点
        gradient_threshold: 算法终止的梯度阈值
        H: 由梯度生成的矩阵, 取代牛顿法计算方向中的hessian矩阵; 
        method: "Newton" & "quasi_Newton_DFP" & "quasi_Newton_BFGS"
        """
        self.x = x0
        self.gradient_threshold = gradient_threshold
        self.method = method
        self.H = H
        
    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 compute_Hessian(self):
        """ 4. 计算当前迭代点的Hessian matrix """
        value = np.mat([[2, 0], [0, 4]], np.float)
        return value
    
    def is_termination(self):
        """ 5. 检验算法是否应该终止 """
        if np.linalg.norm(self.compute_gradient()) < self.gradient_threshold:
            return True
        else:
            return False
        
    def correct_algorithm(self, s, y):
        """ 6. 校正迭代矩阵: BFGS, DFP 
        s: 迭代点位移
        y: 梯度差
        H: 由梯度生成的矩阵, 取代牛顿法计算方向中的hessian矩阵
        """
        if self.method == "quasi_Newton_DFP":
            if s.T * y > 0:
                self.H = self.H - (self.H*y*y.T*self.H)/(y.T*self.H*y) + (s*s.T)/(s.T*y)
        elif self.method == "quasi_Newton_BFGS":
            if y.T*s > 0:
                self.H = self.H - (self.H*s*s.T*self.H)/(s.T*self.H*s) + (y*y.T)/(y.T*s)
        else:
            raise(Exception("No method founded."))
    
    def compute_direction(self):
        """ 7. 计算搜索方向 
        H: 由梯度生成的矩阵, 取代牛顿法计算方向中的hessian矩阵
        """
        if self.method == "Newton":
            value = -np.linalg.inv(self.compute_Hessian()) * self.compute_gradient()
        elif self.method in ["quasi_Newton_DFP", "quasi_Newton_BFGS"]:
            value = -self.H * self.compute_gradient()
        else:
            raise(Exception("No method founded."))
        return value
    
    def compute_step(self, alpha = 1, row = 1/4):
        """ 8. 计算步长 
        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):
        """ 9. 计算下一个迭代点 """
        self.x = self.x + alpha * direction
        
    def draw_result(self, point):
        """ 10. 画出迭代点的变化过程 """
        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_last: 上一个迭代点
        g_last: 上一个迭代点的梯度
        x: 最优解
        """
        iter_point.append(np.linalg.norm(self.x))
        for i in range(max_iter):
            # 算法是否终止
            if self.is_termination():
                break
                
            # 计算方向direction
            direction = self.compute_direction()

            # 计算步长alpha
            alpha = self.compute_step()
            
            # 计算下一个迭代点x
            x_last = self.x.copy()
            g_last = self.compute_gradient().copy()
            self.compute_next_point(alpha, direction)
            
            # 校正迭代矩阵
            if self.method in ["quasi_Newton_BFGS", "quasi_Newton_DFP"]:
                self.correct_algorithm(self.x - x_last, self.compute_gradient() - g_last)
            iter_point.append(np.linalg.norm(self.x))
        return iter_point
        
    
if __name__ == "__main__":
    # init
    params = {"x0": np.mat([[1],[1]]), "gradient_threshold": 1e-6, 
              "method": "quasi_Newton_DFP", "H": np.mat([[1,0],[0,1]], np.float)}
    Newton = Newton_method(**params)
    
    # computing
    iter_point = []
    iter_point = optimal_point = Newton.main(iter_point, max_iter=100)
    
    # draw
    Newton.draw_result(iter_point).show()