数值优化方法实验报告

matlab 计算

共轭梯度法及改进

共轭梯度法

题目

$$\min _{x\in R^5} f(x)=\frac{1}{2}x^TGx-b^Tx,$$
其中 $$G = \begin{pmatrix} 4&1\\ 1&4^2&1\\ &1&4^3&1\\ &&1&4^4&1\\ &&&1&4^5&1\\ &&&&1&4^6&1\\ &&&&&1&4^7&1\\ &&&&&&1&4^8&1\\ &&&&&&&1&4^9&1\\ &&&&&&&&1&4^{10}&1 \end{pmatrix},b=\begin{pmatrix}5\\18\\66\\258\\1026\\4098\\16386\\65538\\262146\\1048577\end{pmatrix}$$

代码

使用共轭梯度法求解该题，代码如下：
main 函数

function GG_main()
G = zeros(10);
for t = 1:10
    G(t,t) = 4^t;
    if t>=2
        G(t,t-1) = 1;
    end
    if t<=9
        G(t,t+1) = 1;
    end
end
disp(G);
b =[5 18 66 258 1026 4098 16386 65538 262146 1048577];
x0 = zeros(1,10);
max_iter = 1000;
fprintf('\n')
fprintf('Conjugate Gradient Methord:\n');
fprintf('=============================');
[y,iter] = frcg(G, b, x0, max_iter );
fprintf('\n')
fprintf(' Iterative number: \n      %d \n',iter);
fprintf(' Solution:\n');
fprintf('   %10.6f',y);
fprintf('\n\n==========================\n')

cg 函数

function [x iter] = cg(G, b, x0, max_iter)
x = x0';
b=b';
tolerance = 1.0e-6;
fprintf('\n x0=');
fprintf('   %10.6f',x0);
%save sb.mat G b x
size(x);
size(b);
size(G);
r = G*x-b;
d=-r;
for k=1:max_iter
    if norm(r,2) <= tolerance
        iter = k-1;
        fprintf('\n Algorithm finds a solution!');
        return
    end
    alpha = (r'*r)/(d'*G*d);
    xx = x+alpha*d;
    rr = r+alpha*G*d;
    beta = (rr'*rr)/(r'*r);
    d = -rr+beta*d;
    x = xx;
    r = rr;
    fprintf('\n  x%d = ',k);
    fprintf('    %10.6f',x);
end
iter = max_iter;
return 

结果如图所示：

改进共轭梯度法

题目

$$\min _{x\in R^5} f(x)=\frac{1}{2}x^TGx-b^Tx,$$
其中 $$G = \begin{pmatrix} 4&1\\ 1&4^2&1\\ &1&4^3&1\\ &&1&4^4&1\\ &&&1&4^5&1\\ &&&&1&4^6&1\\ &&&&&1&4^7&1\\ &&&&&&1&4^8&1\\ &&&&&&&1&4^9&1\\ &&&&&&&&1&4^{10}&1 \end{pmatrix},b=\begin{pmatrix}5\\18\\66\\258\\1026\\4098\\16386\\65538\\262146\\1048577\end{pmatrix}$$ 参考《一种新线性搜索下的FR共轭梯度法》 利用wolfe准则来寻找共轭梯度法的步长 $\alpha$，然后进行迭代，最后获得求解的结果。

代码

main 函数：

function GG_main()
G = zeros(10);
for t = 1:10
    G(t,t) = 4^t;
    if t>=2
        G(t,t-1) = 1;
    end
    if t<=9
        G(t,t+1) = 1;
    end
end
disp(G);
b =[5 18 66 258 1026 4098 16386 65538 262146 1048577];
x0 = zeros(1,10);
max_iter = 1000;
fprintf('\n')
fprintf('Conjugate Gradient Methord:\n');
fprintf('=============================');
[y,iter] = frcg(G, b, x0, max_iter );
fprintf('\n')
fprintf(' Iterative number: \n      %d \n',iter);
fprintf(' Solution:\n');
fprintf('   %10.6f',y);
fprintf('\n\n==========================\n')

frcg 函数：

function [x iter] = frcg(G, b, x0, max_iter)
x = x0';
b=b';
tolerance = 1.0e-6;
fprintf('\n x0=');
fprintf('   %10.6f',x0);
r = G*x-b;
d=-r;
rou = 0.1;
deta = 0.75;
alphamax = 1;
alpha1 = 0;
alpha2 = alphamax;
for k=1:max_iter
    if norm(r,2) <= tolerance
        iter = k-1;
        fprintf('\n Algorithm finds a solution!');
        return
    end
    %temp1 = main(-r'*r,norm(r,2));
    %fprintf('%d',k);
    if k >= 2
        fai1 = 0.5*x'*G*x - b'*x;
        fai1d = -(G*x-b)'*d;
        alpha = 0.618*(alpha2- alpha1)+alpha1;
        fai = 0.5*(x+ alpha*d)'*G*(x+ alpha*d) - b'*(x+ alpha*d);
        while true
            if (fai-fai1) > -rou*alpha*min(fai1d'*d, norm((G*x - b)'*d,2));
                alphab = alpha1 + 0.5*(alpha -alpha1)/(1+(fai1-fai)/((alpha-alpha1)*fai1d));
                alpha2 = alpha;
                alpha = alphab;
            end
            faid =-(G'*(x+alpha*d)-b)'*d;
            if abs(faid'*d) >= deta*min(faid,norm(fai1d,2))
                alphab = alpha + (alpha - alpha1)*faid/(fai1d-faid);
                alpha1 = alpha;
                fai1d = faid;
                fai1 = fai;
                alpha = alphab;
            else                
                alpha = abs(alpha);                
                break
            end
        end
    else
        %fprintf('wo zei')
        alpha = (r'*r)/(d'*G*d);
    end
    fprintf('\n Wolfe alpha: %d',alpha)    
    alpha = (r'*r)/(d'*G*d); 
    fprintf('\n cg alpha: %d',alpha)
    xx = x+alpha*d;
    rr = r+alpha*G*d;
    beta = (rr'*rr)/(r'*r);
    d = -rr+beta*d;
    x = xx;
    r = rr;
    fprintf('\n  x%d = ',k);
    fprintf('    %10.6f',x);
    disp(alpha)
end
end

结果

可以看出wolfe 准则求解获得的 $\alpha$,与共轭梯度法精确线性搜索获得 的步长大小几乎相同。

总结