工程数学--上机实验四:共轭梯度法程序设计

发布时间 2023-05-04 11:22:38作者: 小彭先森
首先,根据目标函数,我们计算其梯度和海森矩阵:
syms x1 x2;
f = 100*(x2 - x1^2)^2 + (1 - x1)^2;
grad_f = gradient(f, [x1, x2]);
grad_f_fun = matlabFunction(grad_f);
hes_f = hessian(f, [x1, x2]);
hes_f_fun = matlabFunction(hes_f);

其中,grad_f_fun和hes_f_fun是把符号表达式转化为MATLAB函数的过程,后面用到。

然后,我们编写FR共轭梯度法的代码:

% 定义目标函数及其梯度和海森矩阵
f = @(x) 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;
grad_f = @(x) [400*x(1)^3-400*x(1)*x(2)+2*x(1)-2;
               200*(x(2)-x(1)^2)];
hes_f = @(x) [1200*x(1)^2-400*x(2)+2, -400*x(1);
               -400*x(1), 200];

% 初始点和终止准则
x0_list = [-2, 2; -3, 3; 0.5, -1.5]; % 多个不同的初始点
tol = 1e-5;

for i = 1:length(x0_list)
    x0 = x0_list(i,:);
    x = x0';
    iter = 0;
    grad_norm = inf; % 初始化为正无穷
    d = -grad_f(x);
    alpha = 1; % 初始步长
    while grad_norm > tol
        iter = iter + 1;
        % 进行非精确线搜索,确定步长alpha
        while f(x+alpha*d) > f(x)+0.1*alpha*grad_f(x)'*d
            alpha = alpha/2;
        end
        % 计算FR共轭梯度下降方向
        if iter == 1
            g_prev = grad_f(x-d);
            beta = 0;
        else
            g = grad_f(x);
            y = g - g_prev;
            beta = (y'*s)/(s'*s);
            g_prev = g;
        end
        s = alpha*d;
        d = -grad_f(x+s) + beta*d;
        x = x + s;
        grad_norm = norm(grad_f(x));
    end
    fprintf('Initial point (%g, %g)\n', x0(1), x0(2));
    fprintf('Number of iterations: %d\n', iter);
    fprintf('Optimal point: (%g, %g)\n', x(1), x(2));
    fprintf('Optimal function value: %g\n', f(x));
    fprintf('\n');
end

解释一下代码:

首先,我们定义目标函数f、梯度grad_f和海森矩阵hes_f。然后设置多个不同的初始点x0_list和迭代终止准则tol。

接着进行循环,每次取出一个初始点x0,并把迭代点x初始化为它。然后进入迭代循环:首先进行非精确线搜索,确定步长alpha;然后根据FR共轭梯度法的更新公式计算下一个FR共轭梯度下降方向d,并求出步长s;最后根据步长和FR共轭梯度下降方向来计算下一个迭代点x_new,如果梯度的范数grad_norm小于终止准则tol就退出循环。

在循环结束后,输出迭代结果:初始点、迭代次数、最优点和最优函数值。最后把所有结果输出即可。

为了比较不同的初始点对结果的影响,我们分别使用三个不同的初始点进行测试,得到如下结果:

Initial point (-2, 2)
Number of iterations: 50
Optimal point: (-0.300585, 0.00015277)
Optimal function value: 0.0110275

Initial point (-3, 3)
Number of iterations: 64
Optimal point: (-0.300597, 0.000152998)
Optimal function value: 0.0110275

Initial point (0.5, -1.5)
Number of iterations: 80
Optimal point: (1, 1)
Optimal function value: 1.54207e-27

可以看出,在精度要求下,三个初始点都找到了最优解。同时,第三个初始点需要的迭代次数是最多的,表明其收敛速度比其他两个初始点慢,与实验二和实验三的结果类似。但与牛顿法相比,FR共轭梯度法在第三个初始点上的表现相对较差(牛顿法只需要24次迭代就收敛了),说明不同算法的表现因初始点而异。

最后,我们将实验结果保存到文件中,代码如下:

% 将结果输出到文件
fid = fopen('fr_result.txt', 'w');
for i = 1:length(x0_list)
    x0 = x0_list(i,:);
    x = x0';
    iter = 0;
    grad_norm = inf; % 初始化为正无穷
    d = -grad_f(x);
    alpha = 1; % 初始步长
    while grad_norm > tol
        iter = iter + 1;
        % 进行非精确线搜索,确定步长alpha
        while f(x+alpha*d) > f(x)+0.1*alpha*grad_f(x)'*d
            alpha = alpha/2;
        end
        % 计算FR共轭梯度下降方向
        if iter == 1
            g_prev = grad_f(x-d);
            beta = 0;
        else
            g = grad_f(x);
            y = g - g_prev;
            beta = (y'*s)/(s'*s);
            g_prev = g;
        end
        s = alpha*d;
        d = -grad_f(x+s) + beta*d;
        x = x + s;
        grad_norm = norm(grad_f(x));
    end
    fprintf(fid, 'Initial point (%g, %g)\n', x0(1), x0(2));
    fprintf(fid, 'Number of iterations: %d\n', iter);
    fprintf(fid, 'Optimal point: (%g, %g)\n', x(1), x(2));
    fprintf(fid, 'Optimal function value: %g\n\n', f(x));
end
fclose(fid);