共轭梯度算法 matlab 实现

function [ fv, bestx, iter_num ] = ConjungateGradient(f, x, x0, epsilon, show_detail ) %conjungate gradient method % Input: % f - syms function % x - row cell arrow for input syms variables % x0 - init point % epsilon - tolerance(公差) % show_detail - a boolean value for whether for print details % Output: % fv - minimum f value % bestx - mimimum point % inter_num - iteration count %% 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) %norm函数 计算范数 %% one-dimensional search %define $x_{k+1} = x_{k} + \lambda d$ xv = xv + lambdas * d; %define $\phi$\phi$ and do 1-dim search phi = subs(f, x, xv); nphi = diff(phi); lambda = solve(nphi); %get rid of complex and minus solution lambda = double(lambda); if length(lambda) > 1 lambda = lambda(abs(imag(lambda)) < 1e-5); %imag 求复数虚数部分 lambda = lambda(lambda > 0); lambda = lambda(1); end % if $\lambdas$ is too samll, stop iteration if lambda < 1e-5 break; end %% update % update $x_{k+1} = x_{k} + \lambds d$ xv = subs(xv, lambdas, lambda); %convert sym to double xv = double(xv); %compute the differential nfv = subs(nf, x, xv); %update counters count = count + 1; k = k + 1; %compute alpha based on FR formula alpha = sumsqr(nfv) / sumsqr(nfv_pre); % show iteration info if show_detail fprintf('Iteration: %d\n', count); fprintf('x(%d) = %s, lambda = %f\n', count, num2str(xv), lambda); fprintf('nf(x) = %s, norm(nf) = %f\n', num2str(double(nfv)), norm(double(nfv))); fprintf('d = %s, alpha = %f\n', num2str(double(d)), double(alpha)); fprintf('\n'); end %update conjungate direction d = -nfv + alpha * d; % save the previous $$\nabla f(x_k)$$ nfv_pre = nfv; % reset the conjungate direction and k if k >= n if k >= 0 k = 0; d = -nfv; end end % while %% output fv = double(subs(f,x,xv)); bestx = double(xv); iter_num = count; end