实验的题目和要求
一、所属课程名称：
最优化方法
二、实验日期：
2010年5月10日~2010年5月15日
三、实验目的
掌握最速下降法，牛顿法和共轭梯度法的算法思想，并能上机编程实现相应的算法。
二、实验要求
用MATLAB实现最速下降法，牛顿法和共轭梯度法求解实例。
四、实验原理
最速下降法是以负梯度方向最为下降方向的极小化算法，相邻两次的搜索方向是互相直交的。牛顿法是利用目标函数)(x
f在迭代点
x处的Taylor展开式作为模型函数，并利用这个二次模型函数的极k
小点序列去逼近目标函数的极小点。共轭梯度法它的每一个搜索方向是互相共轭的，而这些搜索方向
d仅仅是负梯度方向k g-与上一次接
k
待的搜索方向
d的组合。
k
-
1
五．运行及结果如下:
最速下降法：
题目：f=(x-2)^2+(y-4)^2
M文件：
function [R,n]=steel(x0,y0,eps)
syms x;
syms y;
f=(x-2)^2+(y-4)^2;
v=[x,y];
j=jacobian(f,v);
T=[subs(j(1),x,x0),subs(j(2),y,y0)];
temp=sqrt((T(1))^2+(T(2))^2);
x1=x0;y1=y0;
n=0;
syms kk;
while (temp>eps)
d=-T;

最优化算法与matlab应用3：最速下降法
最速下降法是一种沿着N维目标函数的负梯度方向搜索最小值的方法。

（1）算法原理

函数的负梯度表示如下：



搜索步长可调整ak，通常记为 （第k次迭代的步长）。该算法利用一维的线性搜索方法，如二次逼近法，沿着负梯度方向不断搜索函数的较小值，从而找出最优解。

（2）算法步骤



（3）算法实现
%用最速下降法求最优化解
f3= inline('x(1)*(x(1)-5-x(2))+x(2)*(x(2)-4)','x');%目标函数
grad=inline('[2*x(1)-5-x(2),-x(1)+2*x(2)-4]','x'); %目标函数的梯度函数
x0 = [1 4];
TolX = 1e-4; 
TolFun = 1e-9;
MaxIter = 100;
dist0=1;
[xo,fo] = Opt_Steepest(f3,grad,x0,TolX,TolFun,dist0,MaxIter)
运行结果：xo=4.666   4.333    fo=-20.333

function [xo,fo] = Opt_Steepest(f,grad,x0,TolX,TolFun,dist0,MaxIter)
% 用最速下降法求最优化解
%输入：f为函数名 grad为梯度函数
%x0为解的初值 TolX,TolFun分别为变量和函数的误差阈值
%dist0为初始步长 MaxIter为最大迭代次数
%输出: xo为取最小值的点 fo为最小的函数值
% f0 = f(x(0))

%%%%%%判断输入的变量数，设定一些变量为默认值
if nargin < 7
    MaxIter = 100; %最大迭代次数默认为100
end
if nargin < 6
    dist0 = 10; %初始步长默认为10
end
if nargin < 5
    TolFun = 1e-8; %函数值误差为1e-8
end
if nargin < 4
    TolX = 1e-6; %自变量距离误差
end
%%%%%第一步，求解的初值的函数值
x = x0;
fx0 = feval(f,x0);
fx = fx0;
dist = dist0;
kmax1 = 25; %线性搜索法确定步长的最大搜索次数
warning = 0; 
%%%%%迭代计算求最优解

for k = 1: MaxIter
    g = feval(grad,x);
    g = g/norm(g); %求在x处的梯度方向
    %%线性搜索方法确定步长
    dist = dist*2; %令步长为原步长的二倍
    fx1 = feval(f,x-dist*2*g);
    for k1 = 1:kmax1
        fx2 = fx1;
        fx1 = feval(f,x-dist*g);
        if fx0 > fx1+TolFun & fx1 < fx2 - TolFun %fx0 > fx1 < fx2，
            den = 4*fx1 - 2*fx0 - 2*fx2;num = den - fx0 + fx2;  %二次逼近法
            dist = dist*num/den;
            x = x - dist*g; fx = feval(f,x); %确定下一点
            break;
        else
            dist = dist/2;
        end
    end
    if k1 >= kmax1
        warning = warning + 1; %无法确定最优步长
    else
        warning = 0;
    end
    if warning >= 2|(norm(x - x0) < TolX&abs(fx - fx0) < TolFun)
        break;
    end
    x0 = x;
    fx0 = fx;
end
xo = x; fo = fx;
if k == MaxIter
    fprintf('Just best in %d iterations',MaxIter);
end

题目和要求
最速下降法是以负梯度方向最为下降方向的极小化算法，相邻两次的搜索方向是互相直交的。牛顿法是利用目标函数f(x)在迭代点xk处的Taylor展开式作为模型函数，并利用这个二次模型函数的极
小点序列去逼近目标函数的极小点。共轭梯度法它的每一个搜索方向是互相共轭的，而这些搜索方向dk仅仅是负梯度方向 gk与上一次接待的搜索方向dk 1的组合。
运行及结果如下:
最速下降法：
题目：f=(x-2)^2+(y-4)^2
M文件：
function [R,n]=steel(x0,y0,eps)
syms x;
syms y;
f=(x-2)^2+(y-4)^2;
v=[x,y];
j=jacobian(f,v);
T=[subs(j(1),x,x0),subs(j(2),y,y0)];
temp=sqrt((T(1))^2+(T(2))^2);
x1=x0;y1=y0;
n=0;
syms kk;
while (temp>eps)
d=-T;
f1=x1+kk*d(1);f2=y1+kk*d(2);
fT=[subs(j(1),x,f1),subs(j(2),y,f2)];
fun=sqrt((fT(1))^2+(fT(2))^2);
Mini=Gold(fun,0,1,0.00001);
x0=x1+Mini*d(1);y0=y1+Mini*d(2);
T=[subs(j(1),x,x0),subs(j(2),y,y0)];
temp=sqrt((T(1))^2+(T(2))^2);
x1=x0;y1=y0;
n=n+1;
end
R=[x0,y0]
调用黄金分割法：

最近在课程实践作业中，要求用最速下降法、牛顿法和拟牛顿法三种方法求解高维一致凸二次优化问题的极小值，网上看到的大部分程序都是手动求好了凸二次函数 f 的偏导然后带进去计算，这样的话限制死了维数和次数，也让程序显得比较笨拙，因此就自己用python从零实现了一下，由于要的急也还有很多改进空间吧…这种优化问题其实用 matlab 会比较方便，因此在 python 里想的也是借鉴 matlab中的符号计算体系，所以基于 sympy 库去实现的。
最速下降法
理论部分：

代码部分：
import sympy as sy #sympy是Python的科学计算库，引入符号计算体系
import matplotlib.pyplot as plt #导入matplot库，画图所需
import math
import time

def steepest_descent(f, x): #输入函数f，初始点x，根据最速下降法求解无约束最优化问题
    n = len(x) #维数n
    diff_f = [] #存储梯度向量
    sym_x = [] #存储自变量的符号形式
    for i in range(n): 
        sym = 'x' + str(i+1)
        sym_x.append(sy.Symbol(sym)) #sym_x = [x1,x2,……,xn]
        diff = sy.diff(f, sym_x[i]) #计算梯度
        diff_f.append(diff)
    
    error = 1 #误差
    time = 0 #迭代次数
    list_time = [] #存储画图所需横坐标
    list_error = [] #存储画图所需纵坐标
    while(error > 1e-6) and (time < 100) : #设置终止循环条件
        d = [] #搜索方向d(k)
        sub_d = {}
        for i in range(n):
            sub_d[sym_x[i]] = x[i]
        for i in range(n):
            d.append( - diff_f[i].evalf(subs=sub_d)) #相当于把x1等参数换成了具体数值，代入梯度计算
        print("当前搜索方向为：", d)
        
        m = sy.Symbol('m') #步长变量alpha
        func = f
        for i in range(n):
            func = func.subs(sym_x[i], x[i] + m * d[i]) # 实现f(x) = f(x + m * d)，相当于把x1等参数换成了含m变量，此时func是关于m的函数
        best_m = sy.solve(sy.diff(func)) #解得最优步长m(k)
        print("当前最优步长为：", best_m)
        
        x1 = [] #计算下一个点
        for i in range(n):
            x1.append(x[i] + best_m[0] * d[i]) #x(k+1) = x(k) + m(k) * d(k)
        x = x1 #迭代
        
        res = 0 #计算误差
        sub_x = {}
        for i in range(n):
            sub_x[sym_x[i]] = x[i]
        for i in range(n):
            res += diff_f[i].evalf(subs=sub_x) ** 2 #同样是传入具体数值代替参数进行计算
        error = math.sqrt(res)
        print("当前误差：", error)
        list_error.append(error)
        time += 1
        list_time.append(time)
    print("循环次数为：",time)
    plt.plot(list_time, list_error, 'o-', color = 'r') #设置画图参数
    plt.xlabel("Number of iterations") #横坐标名字
    plt.ylabel("The error") #纵坐标名字
    plt.show()
    return x
        
        
t0 = time.clock()
x1, x2, x3, x4, x5, x6, x7, x8, x9, x10 = sy.symbols('x1, x2, x3, x4, x5, x6, x7, x8, x9, x10')
f = x1**2 + 2*x2**2 + 3*x3**2 + 4*x4**2 + 5*x5**2 + 6*x6**2 + 7*x7**2 + 8*x8**2 + 9*x9**2 + 10*x10**2 - x1 - 9*x2 - 8*x3 - 7*x4 - 6*x5 - 2*x6 - 4*x7 - 4*x8 - x9 - 3*x10
x = steepest_descent(f, [100,0,1,0,100,0,1,0,100,0])
print("解为：", x)
print("所用时间为", time.clock()-t0,"秒")

运行结果为：



牛顿法
理论部分：



代码部分：
import sympy as sy #sympy是Python的科学计算库，引入符号计算体系
import matplotlib.pyplot as plt #导入matplot库，画图所需
import numpy as np #导入numpy库，矩阵运算所需
import math
import time
import copy

def newton(f, x): #输入函数f，初始点x，根据牛顿法求解无约束最优化问题
    n = len(x) #维数n
    diff_f = [] #存储梯度向量
    sym_x = [] #存储自变量的符号形式
    for i in range(n): 
        sym = 'x' + str(i+1)
        sym_x.append(sy.Symbol(sym)) #sym_x = [x1,x2,……,xn]
        diff = sy.diff(f, sym_x[i]) #计算梯度
        diff_f.append(diff)

    hesse_matrix = [] #存储海森矩阵
    for i in range(n):
        row = []#当前行的梯度
        for j in range(n):
            row.append(sy.diff(diff_f[i], sym_x[j])) #计算梯度
        hesse_matrix.append(row)       
        
    error = 1 #误差
    time = 0 #迭代次数
    list_time = [0] #存储画图所需横坐标
    list_error = [1] #存储画图所需纵坐标
    while(error > 1e-6) and (time < 100) : #设置终止循环条件
        g = copy.deepcopy(diff_f) #计算梯度向量g(k)
        H = copy.deepcopy(hesse_matrix) #计算海森矩阵H(k)
        sub_g = {}
        for i in range(n):
            sub_g[sym_x[i]] = x[i]
        for i in range(n):
            g[i] = float(diff_f[i].evalf(subs=sub_g))
            for j in range(n):
                H[i][j] = float(hesse_matrix[i][j].evalf(subs=sub_g))
        print("当前g(k)：", g)
        print("当前H(k)：", H)
        
        d = - np.dot(np.linalg.inv(np.array(H)), np.array(g)) #牛顿方向d(k) = -H(k)^-1 * g(k),inv()实现求逆
        x += d #迭代，x(k+1) = x(k) + d
        
        res = 0 #计算误差
        sub_x = {}
        for i in range(n):
            sub_x[sym_x[i]] = x[i]
        for i in range(n):
            res += diff_f[i].evalf(subs=sub_x) ** 2
        error = math.sqrt(res)
        print("当前误差：", error)
        list_error.append(error)
        time += 1
        list_time.append(time)
    print("循环次数为：",time)
    plt.plot(list_time, list_error, 'o-', color = 'r') #设置画图参数
    plt.xlabel("Number of iterations") #横坐标名字
    plt.ylabel("The error") #纵坐标名字
    plt.show()
    return x

拟牛顿法：
理论部分：



代码部分：
import sympy as sy #sympy是Python的科学计算库，引入符号计算体系
import matplotlib.pyplot as plt #导入matplot库，画图所需
import numpy as np #导入numpy库，矩阵运算所需
import math
import time
import copy

def quasi_newton(f, x): #输入函数f，初始点x，根据DFP拟牛顿法求解无约束最优化问题
    n = len(x) #维数n
    diff_f = [] #存储梯度向量
    sym_x = [] #存储自变量的符号形式
    for i in range(n): 
        sym = 'x' + str(i+1)
        sym_x.append(sy.Symbol(sym)) #sym_x = [x1,x2,……,xn]
        diff = sy.diff(f, sym_x[i]) #计算梯度
        diff_f.append(diff)
    
    D = np.zeros((n, n)) 
    for i in range(n):
        D[i][i] = 1 #D(0) = I 单位矩阵
    error = 1 #误差
    time = 0 #迭代次数
    list_time = [] #存储画图所需横坐标
    list_error = [] #存储画图所需纵坐标
    while(error > 1e-6) and (time < 100) : #设置终止循环条件
        g = copy.deepcopy(diff_f) #计算梯度向量g(k)
        sub_g = {}
        for i in range(n):
            sub_g[sym_x[i]] = x[i]
        for i in range(n):
            g[i] = float(diff_f[i].evalf(subs=sub_g))
        
        d = - D@g #d(k) = - D(k) * g(k),@是numpy中另一种矩阵乘法的实现
        print("当前搜索方向为：", d)
        
        m = sy.Symbol('m') #步长变量
        func = f
        for i in range(n):
            func = func.subs(sym_x[i], x[i] + m * d[i]) # f(x) = f(x + m * d)
        best_m = sy.solve(sy.diff(func)) #解得最优步长m(k)
        print("当前最优步长为：", best_m)
        
        s = float(best_m[0]) * d #s(k) = m(k) * d(k)
        x = x + s #x(k+1) = x(k) + s(k)
        
        g1 = copy.deepcopy(g) #计算梯度向量g(k+1)
        sub_g = {}
        for i in range(n):
            sub_g[sym_x[i]] = x[i]
        for i in range(n):
            g1[i] = float(diff_f[i].evalf(subs=sub_g))
        
        y = np.array(g1) - np.array(g) #y(k) = g(k+1) - g(k)
        s.shape, y.shape = (n,1), (n,1) #计算一维数组转置需要补充维度
        # D(k+1) = D(k) + s(k)*s(k)^T / s(k)^T*y(k) - D(k)*y(k)*y(k)^T*D(k) / y(k)^T*D(k)*y(k)
        D += (np.dot(s, s.transpose()) / np.dot(s.transpose(), y)) - (np.dot(np.dot(np.dot(D, y), y.transpose()), D) / np.dot(np.dot(y.transpose(), D), y))
        
        res = 0 #计算误差
        for i in range(n):
            res += g1[i]** 2
        error = math.sqrt(res)
        print("当前误差：", error)
        list_error.append(error)
        time += 1
        list_time.append(time)
    print("循环次数为：",time)
    plt.plot(list_time, list_error, 'o-', color = 'r') #设置画图参数
    plt.xlabel("Number of iterations") #横坐标名字
    plt.ylabel("The error") #纵坐标名字
    plt.show()
    return x

存在的问题
要说略显笨拙的地方首先就是符号计算吧，由于第一次像这样大面积使用sympy库，我总感觉同样的思路有更简单的函数什么的，但又不太熟悉；其次就是函数 f 还是得自己敲，如果维数高一点也很麻烦，但输入系数矩阵 A 我感觉也同样麻烦…
另外，程序中的示例函数是十维二次函数，更高维度自然是没有问题，但是更高次数，比如三次、四次函数在


这里会有点问题，这里输出的best_m是个列表，也可能解为虚数，虚数显然没办法进行下一步计算，会报错。