MATLAB中使用基础优化算法求解二次函数问题

啊第一次在csdn上写东西，在前些时候优化方法的计算实习老师留了一些上机作业，第一道题目是用不同的优化算法结合最速下降法求解二次函数。

二次函数及梯度的生成

先贴代码/matlab:

function Function_value = Quadratic_function(x)
%Quadratic_function(x)函数用来对给定的x返回二次函数值
n = 146;
a = unidrnd(10,n,1);
G = a*a'+unidrnd(2)*eye(n);
b = 0.5*G*ones(n,1);
Function_value = 1/2*(x')*G*x-(b')*x;

我们使用unidrnd函数生成随机整数的列向量，我们的基础二次函数的形式是：

以及写出的梯度函数为之后做准备，不过这个函数的梯度特别好求，写不写都可。

function Gradient_value = Function_gradient(x)
%我们将要计算函数Quadratic_function对146个变量的导数
Gradient_value = G*x - b;

最速下降法加精确线搜索

function  xk =  Steepest_plus_Exact_line_search()
%我们采用精确线搜索加最速下降法
x0 = zeros(146,1);
eps = 0.001;
n = 146;
a = unidrnd(10,n,1);
G = a*a'+unidrnd(2)*eye(n);
b = 0.5*G*ones(n,1);
%做循环的准备条件
k = 0;xk = x0;
res = norm(Function_gradient(xk));
while res > eps && k<=10000
%dk
gk = G*xk-b;
dk = -gk;
%alphak
alphak = (gk'*gk)/(gk'*G*gk);
xk = xk + alphak*dk;
%更新k
k = k+1;
res = norm(gk);
%打印循环情况
fprintf('在%d-th轮，梯度的值为——%f\n',k,res);
end

啊最速下降法的收敛速度是最慢的。

牛顿法加精确线搜索

function [xk,k] = Newtown_plus_Exact_line_search()
%我们采用精确线搜索加最速下降法
x0 = zeros(146,1);
eps = 0.001;
n = 146;
a = unidrnd(10,n,1);
G = a*a'+unidrnd(2)*eye(n);
b = 0.5*G*ones(n,1);
%做循环的准备条件
gk = Function_gradient(x0);
res = norm(gk);
k = 0;
xk = x0;
G_inverse = inv(G);
alphak = 1;
while res > eps && k<=10000
%dk
dk = -G_inverse*gk;
gk = G*xk-b;
alphak = -(dk'*gk)/(dk'*G*dk);
%xk
xk = xk + alphak * dk;
%f = Quadratic_function(xk);
k = k + 1;
res = norm(gk);
%打印循环结果
fprintf('在%d-th轮，梯度的值为——%f\n',k,res);
end

然后牛顿法是最快的，因为实际上对于二次函数问题来说的话，牛顿法可以一步达到最小值点。

共轭梯度法

function [x_star]=Conjugate_gradient_method()
%在该脚本下我们使用共轭梯度法进行求解
%给定x0，eps
x0 = zeros(146,1);
eps = 0.001;
%为方便使用将条件列出：
n = 146;
a = unidrnd(10,n,1);
G = a*a'+unidrnd(2)*eye(n);
b = 0.5*G*ones(n,1);
gradk =  G*x0 - b;
k=0;
%选择初始点，最大迭代次数
xk=x0;
maxit=10000;
while 1
%gradk的二范数
res=norm(gradk);
%第k次迭代
k=k+1;
%算法停止条件
if res<=eps||k>maxit

break;
end
%确定下降方向
dk=-gradk;
%计算步长
alphak=-gradk'*dk/(dk'*G*dk);
%迭代核心
xk=xk+alphak*dk;
%计算Beta
grad_temp=gradk;
gradk=G*xk-b;
Betak=gradk'*gradk/(grad_temp'*grad_temp);
dk=-gradk+Betak*dk;
fprintf('At the %d-th iteration,the residual is  %f\n',k,res);
end
x_star=xk;

共轭梯度达到最小值不超过n步，但迭代一段时间之后，正交性会消失，实际步数可能会比n多。

BFGS加精确线搜索

function[x_star]= BFGS_plus_Exact_line_search()
%在该脚本下我们使用BFGS加精确线搜索求解问题
%给定x0，eps
x0 = zeros(146,1);
eps = 0.001;
xk=x0;
%为方便使用将条件列出：
n = 146;
a = unidrnd(10,n,1);
G = a*a'+unidrnd(2)*eye(n);
b = 0.5*G*ones(n,1);
%gk，hess阵
gradk = Function_gradient(xk);hess = G;
res=norm(gradk);k=0;Hk=eye(146);
%开始循环
while res>eps&&k<=10000
dk=-Hk*gradk;
%exact line search
alphak=-dk'*gradk/(dk'*hess*dk);
%保存xk
xkk=xk;
%计算下一点
xk=xk+alphak*dk;
%两点差
deltak=xk-xkk;
%保存xk的梯度
gradkk=gradk;
%unidrnd生成的数是随机的,要用第一次hess阵和b
gradk=hess*xk-b;
%两点梯度差
yk=gradk-gradkk;
%bfgs公式
Hk=(eye(146)-(deltak*yk')/(deltak'*yk))*Hk*(eye(146)-(yk*deltak')/(deltak'*yk))+(deltak*deltak')/(deltak'*yk);
k=k+1;
%下一点梯度范数
res=norm(gradk);
fprintf('At the %d-th iteration,the residual is:----%f\n',k,res)
end
x_star=xk;
end

最后通过这四种方法计算出来的结果应该均是收敛到0.5，所以到此就做完了这道题！

1. 不过有一说一，在做完之后看代码确实有好多问题，其实我觉得首先梯度文件就不用写，每次要求梯度的时候直接用Gx-b就很方便，而且写的过程中，也要一直用初始的G和b，但是如果多次调用function_value函数的话就每次随机得到的也是不同的矩阵，所以就会出问题，也是自己遇到的问题之一。
2. 第二个问题是写这种算法，很容易思路混乱，所以下次如果再有类似的情况，希望自己可以先将算法一步一步写出来，这样的话编程序的时候也思路比较清楚，而且修改代码也更容易。

啊到这第一篇就写完了！！完事 收工！！

• 点赞 1
• 收藏
• 分享
• 文章举报