共轭梯度法求解无约束最优化问题

一、共轭梯度法基本原理

1. 算法核心思想

共轭梯度法是一种用于求解大规模无约束优化问题的高效迭代算法，特别适用于目标函数为二次型或大规模稀疏系统的情况。其核心思想是通过构造一组共轭方向，在每一步迭代中沿这些方向进行精确线搜索，从而在有限步内收敛到最优解。

2. 数学基础

对于二次型目标函数：

\(f(x)=\frac{1}{2}x^TAx+b^Tx+c\)

其中 A是对称正定矩阵

共轭方向定义：

\(d_i^TAd_j=0,i\neq{j}\)

3. 算法优势

方法	存储需求	计算复杂度	收敛速度	适用问题规模
梯度下降法	O(n)	O(n)	线性收敛	中小规模
牛顿法	O(n²)	O(n³)	二次收敛	小规模
共轭梯度法	O(n)	O(n²)	超线性收敛	大规模

二、算法推导与实现

1. 标准共轭梯度算法（Fletcher-Reeves）

function [x_opt, f_opt, iter] = conjugateGradient(fun, grad, x0, tol, max_iter)
    % 输入参数:
    %   fun: 目标函数句柄
    %   grad: 梯度函数句柄
    %   x0: 初始点
    %   tol: 收敛容差
    %   max_iter: 最大迭代次数
    
    x = x0;
    g = grad(x);          % 初始梯度
    d = -g;               % 初始搜索方向
    f = fun(x);
    
    for iter = 1:max_iter
        % 1. 线搜索确定步长α
        alpha = lineSearch(fun, grad, x, d);
        
        % 2. 更新迭代点
        x_new = x + alpha * d;
        g_new = grad(x_new);
        f_new = fun(x_new);
        
        % 3. 检查收敛条件
        if norm(g_new) < tol
            x_opt = x_new;
            f_opt = f_new;
            return;
        end
        
        % 4. 计算共轭系数β (Fletcher-Reeves公式)
        beta = (g_new' * g_new) / (g' * g);
        
        % 5. 更新搜索方向
        d = -g_new + beta * d;
        
        % 6. 更新变量
        x = x_new;
        g = g_new;
        f = f_new;
    end
    
    % 达到最大迭代次数仍未收敛
    warning('达到最大迭代次数，可能未收敛');
    x_opt = x;
    f_opt = f;
    iter = max_iter;
end

2. 精确线搜索实现

function alpha = exactLineSearch(fun, grad, x, d)
    % 对于二次函数 f(x) = 1/2*x'*A*x + b'*x
    % α = argmin f(x + αd) = - (g'*d) / (d'*A*d)
    
    g = grad(x);
    A = hessianApproximation(fun, x);  % 近似Hessian矩阵
    
    numerator = g' * d;
    denominator = d' * A * d;
    
    if denominator <= eps
        alpha = 0;  % 避免除以零
    else
        alpha = -numerator / denominator;
    end
end

3. 非精确线搜索（Wolfe条件）

function alpha = wolfeLineSearch(fun, grad, x, d, alpha_max)
    % Wolfe条件线搜索
    c1 = 1e-4;  % Armijo条件常数
    c2 = 0.9;   % 曲率条件常数
    
    alpha = alpha_max;
    phi0 = fun(x);
    dphi0 = grad(x)' * d;
    
    while true
        x_new = x + alpha * d;
        phi = fun(x_new);
        dphi = grad(x_new)' * d;
        
        % Armijo条件
        armijo = phi <= phi0 + c1 * alpha * dphi0;
        
        % 曲率条件
        curvature = dphi >= c2 * dphi0;
        
        if armijo && curvature
            break;
        elseif ~armijo
            alpha = alpha * 0.5;  % 减小步长
        else
            alpha = alpha * 1.5;  % 增大步长
        end
        
        if alpha < 1e-10  % 防止过小的步长
            alpha = 0;
            break;
        end
    end
end

三、算法变体与改进

1. Polak-Ribière-Polyak (PRP) 公式

% 替代Fletcher-Reeves的β计算公式
beta_PRP = (g_new' * (g_new - g)) / (g' * g);
d = -g_new + beta_PRP * d;

2. 重启策略（Restart Strategy）

% 每n步重启一次（n为变量维度）
if mod(iter, length(x0)) == 0
    d = -g_new;  % 重置为负梯度方向
end

3. 预处理共轭梯度法

function d = preconditionedDirection(g, M_inv)
    % 预处理技术：M为对称正定预处理矩阵
    z = M_inv * g;  % 近似解方程 Mz = g
    beta = (z' * z) / (d_prev' * d_prev);  % 类似FR公式
    d = -z + beta * d_prev;
end

四、MATLAB实现示例

1. 求解Rosenbrock函数最小值

% 定义Rosenbrock函数及其梯度
rosenbrock = @(x) (1 - x(1))^2 + 100*(x(2) - x(1)^2)^2;
rosenbrock_grad = @(x) [-2*(1 - x(1)) - 400*x(1)*(x(2) - x(1)^2); 
                        200*(x(2) - x(1)^2)];

% 初始点和参数设置
x0 = [-1.5; 2.0];
tol = 1e-6;
max_iter = 1000;

% 调用共轭梯度法
[x_opt, f_opt, iter] = conjugateGradient(...
    rosenbrock, rosenbrock_grad, x0, tol, max_iter);

% 显示结果
fprintf('最优解: x1 = %.6f, x2 = %.6f\n', x_opt(1), x_opt(2));
fprintf('最小值: f(x) = %.6f\n', f_opt);
fprintf('迭代次数: %d\n', iter);

2. 求解线性方程组（Ax=b）

% 生成对称正定矩阵A和右端项b
n = 100;  % 问题规模
A = gallery('poisson', sqrt(n));  % Poisson方程离散矩阵
b = randn(n, 1);

% 定义二次型函数和梯度
quad_fun = @(x) 0.5*x'*A*x - b'*x;
quad_grad = @(x) A*x - b;

% 初始点
x0 = zeros(n, 1);

% 调用共轭梯度法
tic;
[x_sol, f_min, iter] = conjugateGradient(...
    quad_fun, quad_grad, x0, 1e-8, n);
toc;

% 验证结果
fprintf('残差范数: %.2e\n', norm(A*x_sol - b));
fprintf('迭代次数: %d (理论最大值=%d)\n', iter, n);

五、性能分析与优化

1. 收敛性分析

对于二次型函数 \(f(x)=\frac{1}{2}x^TAx−b^Tx\)：

• 精确算术运算下，最多n步收敛（n为变量维度）
• 实际计算中，由于舍入误差，通常需要重启策略

收敛条件数依赖：

\(κ(A)=\frac{λmax}{λmin}\)

收敛速度与条件数密切相关

2. 预处理技术

预处理方法	适用场景	效果
对角预处理	变量尺度差异大	改善条件数
不完全Cholesky	稀疏矩阵	保持稀疏性
几何多重网格	PDE相关问题	大幅加速收敛

3. 并行化实现

% 并行计算梯度
parfor i = 1:n
    grad(i) = computePartialGradient(x, i);
end

% 分布式存储方向向量
spmd
    local_d = getLocalDirection();
end

六、工程应用案例

1. 有限元结构分析

% 求解弹性力学方程 K*u = F
K = assembleStiffnessMatrix(mesh);  % 组装刚度矩阵
F = assembleLoadVector(mesh);      % 组装载荷向量

% 使用预处理共轭梯度法
u0 = zeros(size(F));
u = pcgSolver(K, F, u0, struct('M', diag(diag(K))));  % 对角预处理

2. 机器学习模型训练

% 逻辑回归参数优化
X = load('data.csv');  % 特征矩阵
y = load('labels.csv'); % 标签向量

% 正则化逻辑损失函数
fun = @(w) logisticLoss(X, y, w, lambda);
grad = @(w) logisticGrad(X, y, w, lambda);

% 使用共轭梯度法训练
w0 = zeros(size(X, 2), 1);
w_opt = conjugateGradient(fun, grad, w0, 1e-6, 500);

3. 电磁场计算

% 求解Maxwell方程组
A = assembleMaxwellMatrix(mesh, omega);  % 复对称矩阵
E0 = initialFieldGuess();

% 使用复共轭梯度法
E = complexCG(A, E0, struct('tol', 1e-8));

参考代码 共轭梯度法求解无约束最优化问题 www.youwenfan.com/contentcnn/83629.html

七、高级主题扩展

1. 非线性共轭梯度法

function beta = daiYuanFormula(g_new, g_old, d_old)
    % Dai-Yuan公式（适用于非二次函数）
    numerator = g_new' * (g_new - g_old);
    denominator = d_old' * (g_new - g_old);
    beta = numerator / denominator;
end

2. 分布式共轭梯度法

% 使用MPI并行计算
MPI_Init();
comm = MPI_COMM_WORLD;
rank = MPI_Comm_rank(comm);

% 各进程计算局部梯度
local_grad = computeLocalGradient(local_x);

% 全局归约求和
global_grad = MPI_Allreduce(local_grad, comm, 'SUM');

% 广播新方向
MPI_Bcast(d, comm);

3. GPU加速实现

% 使用GPU数组加速计算
x_gpu = gpuArray(x0);
g_gpu = gpuArray(zeros(size(x0)));

% 在GPU上计算梯度
g_gpu = arrayfun(@computeGradient, x_gpu);

% 使用cuBLAS进行矩阵运算
alpha = cublasDotProduct(g_gpu, d_gpu) / ...
        cublasDotProduct(d_gpu, A*d_gpu);