Matlab实现FR共轭梯度法 |
您所在的位置：网站首页 › matlab共轭梯度法函数 › Matlab实现FR共轭梯度法 |
Matlab实现FR共轭梯度法
阅读数 2,436
 +5 前一段时间学习了无约束最优化方法，今天用Matlab实现了求解无约束最优化问题的FR共轭梯度法。关于共轭梯度法的理论介绍，请参考我的另一篇文章 无约束最优化方法学习笔记。 文件testConjungateGradient.m用于测试共轭梯度法函数。测试文件需要定义函数$f$和自变量$x$，给定迭代初值$x_0$和允许误差$\epsilon$。函数设置了show_detail变量用于控制是否显示每一步的迭代信息。 % test conjungate gradient method % by TomHeaven, [email protected], 2015.08.25 %% define function and variable syms x1 x2; %f = xs^2+2*ys^2-2*xs*ys + 2*ys + 2; f = (x1-1)^4 + (x1 - x2)^2; %f = (1-x1)^2 + 2*(x2 - x1^2)^2; x = {x1, x2}; % initial value x0 = [0 0]; % tolerance epsilon = 1e-1; %% call conjungate gradient method show_detail = true; [bestf, bestx, count] = conjungate_gradient(f, x, x0, epsilon, show_detail); % print result fprintf('bestx = %s, bestf = %f, count = %d\n', num2str(bestx), bestf, count);文件conjungate_gradient.m是共轭梯度法的实现函数。变量nf表示函数$f$的梯度$\nabla f$（梯度的希腊字母是nabla，故用nf）。 function [fv, bestx, iter_num] = conjungate_gradient(f, x, x0, epsilon, show_detail) %% conjungate gradient method % by TomHeaven, [email protected], 2015.08.25 % Input: % f - syms function % x - row cell arrow for input syms variables % $x_0$ - init point % epsilon - tolerance % show_detail - a boolean value for wether to print details % Output: % fv - minimum f value % bestx - mimimum point % iter_num - iteration count %% init syms lambdas % suffix s indicates this is a symbol variable % n is the dimension n = length(x); % compute differential of function f stored in cell nf nf = cell(1, n); % using row cells, column cells will result in error for i = 1 : n nf{i} = diff(f, x{i}); end % $\nabla f(x_0)$ nfv = subs(nf, x, x0); % init $\nabla f(x_k)$ nfv_pre = nfv; % init count, k and xv for x value. count = 0; k = 0; xv = x0; % initial search direction d = - nfv; % show initial info if show_detail fprintf('Initial:\n'); fprintf('f = %s, x0 = %s, epsilon = %f\n\n', char(f), num2str(x0), epsilon); end %% loop while (norm(nfv) > epsilon) %% one-dimensional search % define $x_{k+1} = x_{k} + \lambda d$ xv = xv+lambdas*d; % define $\phi$ and do 1-dim search phi = subs(f, x, xv); nphi = diff(phi); % $\nabla \phi$ lambda = solve(nphi); % get rid of complex and minus solution lambda = double(lambda); if length(lambda) > 1 lambda = lambda(abs(imag(lambda)) 0); lambda = lambda(1); end % if $\lambda$ is too small, stop iteration if lambda = n if k >= n k = 0; d = - nfv; end end % while %% output fv = double(subs(f, x, xv)); bestx = double(xv); iter_num = count; end运行testConjungateGradient后输出结果如下： >> testConjungateGradient Initial: f = (x1 - x2)^2 + (x1 - 1)^4, x0 = 0 0, epsilon = 0.100000 Iteration: 1 x(1) = 0.41025 0, lambda = 0.102561 nf(x) = 1.08e-16 -0.82049, norm(nf) = 0.820491 d = 4 0, alpha = 0.042075 Iteration: 2 x(2) = 0.52994 0.58355, lambda = 0.711218 nf(x) = -0.52265 0.10721, norm(nf) = 0.533528 d = 0.1683 0.82049, alpha = 0.422831 Iteration: 3 x(3) = 0.63914 0.56115, lambda = 0.208923 nf(x) = -0.031994 -0.15597, norm(nf) = 0.159223 d = 0.52265 -0.10721, alpha = 0.089062 Iteration: 4 x(4) = 0.76439 0.79465, lambda = 1.594673 nf(x) = -0.11285 0.060533, norm(nf) = 0.128062 d = 0.078542 0.14643, alpha = 0.646892 Iteration: 5 x(5) = 0.79174 0.77998, lambda = 0.242379 nf(x) = -0.012614 -0.023517, norm(nf) = 0.026686 d = 0.11285 -0.060533, alpha = 0.043425 bestx = 0.79174 0.77998, bestf = 0.002019, count = 5修改允许误差为 epsilon = 1e-8;则可以得到更加精确的结果： Iteration: 6 x(6) = 0.9026 0.9122, lambda = 6.329707 nf(x) = -0.022884 0.019188, norm(nf) = 0.029864 d = 0.017515 0.020888, alpha = 1.252319 Iteration: 7 x(7) = 0.90828 0.90744, lambda = 0.247992 nf(x) = -0.0014077 -0.0016788, norm(nf) = 0.002191 d = 0.022884 -0.019188, alpha = 0.005382 Iteration: 8 x(8) = 0.97476 0.97586, lambda = 43.429293 nf(x) = -0.0022668 0.0022025, norm(nf) = 0.003161 d = 0.0015309 0.0015756, alpha = 2.080989 Iteration: 9 x(9) = 0.97533 0.97531, lambda = 0.249812 nf(x) = -2.9597e-05 -3.0461e-05, norm(nf) = 0.000042 d = 0.0022668 -0.0022025, alpha = 0.000181 Iteration: 10 x(10) = 0.99709 0.99712, lambda = 725.188481 nf(x) = -5.2106e-05 5.2008e-05, norm(nf) = 0.000074 d = 3.0006e-05 3.0063e-05, alpha = 3.004594 Iteration: 11 x(11) = 0.9971 0.9971, lambda = 0.249997 nf(x) = -4.8571e-08 -4.8663e-08, norm(nf) = 0.000000 d = 5.2106e-05 -5.2008e-05, alpha = 0.000001 Iteration: 12 x(12) = 0.99992 0.99992, lambda = 57856.826721 nf(x) = -9.3751e-08 9.3748e-08, norm(nf) = 0.000000 d = 4.8616e-08 4.8617e-08, alpha = 3.718503 Iteration: 13 x(13) = 0.99992 0.99992, lambda = 0.250000 nf(x) = -1.1858e-12 -1.1855e-12, norm(nf) = 0.000000 d = 9.3751e-08 -9.3748e-08, alpha = 0.000000 bestx = 0.99992 0.99992, bestf = 0.000000, count = 13这与问题的最优解$(1,1)^T$已经非常接近了。 算法实现没有经过大量测试，实际使用可能会有BUG。这里只是用于说明基本实现原理，有兴趣的读者可以在此基础上改进。 ×用微信扫描并分享 |
【本文地址】
公司简介
联系我们
今日新闻 |
点击排行 |
|
推荐新闻 |
图片新闻 |
|
专题文章 |