截斷共軛梯度法

考慮信賴域子問題： 其中 是目標函數，$\nabla f(x), \nabla^2 f(x)$ 表示 的梯度與海瑟矩陣。注意，當 時，信賴域子問題就等同于求解牛頓方程。

這里，實現截斷共軛梯度法 (Steihaug-Toint Conjugate gradient, ST-CG 方法)來求解上述信賴域子問題。

當約束不存在時(即 )，共軛梯度法通過求解一系列共軛方向可以快速求解對應的無約束二次優化問題。具體地，對于問題 ，給定初始 , , ，共軛梯度法的迭代格式為：

當信賴域約束存在時，共軛梯度法得到的解不能保證落在可行域內。因此，截斷共軛梯度法則在共軛梯度法增加兩條額外的終止條件，用于處理負曲率 或超過信賴域半徑 的情況。

目錄

初始化和迭代準備

輸入信息：迭代點 ，梯度 grad，海瑟矩陣 hess，信賴域半徑 ， 包含算法參數的結構體 opts 。 輸出信息：信賴域子問題的解 (也就是迭代點 處的下降方向)， Heta 為海瑟矩陣在點 作用在方向 上的結果，即 ，記錄迭代信息的結構體 out 和記錄退出原因的標記 stop_tCG 。function [eta, Heta, out, stop_tCG] ...

= tCG(x, grad, hess, Delta, opts)

從輸入的結構體 opts 中讀取參數或采取默認參數。

opts.kappa：牛頓方程求解精度的參數

theta：牛頓方程求精精度的參數

opts.maxit：最大迭代次數，對于共軛梯度法而言默認等于變量的維度

opts.minit：最小迭代次數if ~isfield(opts,'kappa'); opts.kappa = 0.1; end

if ~isfield(opts,'theta'); opts.theta = 1; end

if ~isfield(opts,'maxit'); opts.maxit = length(x); end

if ~isfield(opts,'minit'); opts.minit = 5; end

% 參數復制。

theta = opts.theta;

kappa = opts.kappa;

初始化， 為優化變量，初始為全 0向量：$\mathbf{0}$。 初始化為目標函數的梯度 。eta = zeros(size(x));

Heta = eta;

r = grad;

e_Pe = 0;

r_r = r'*r;

norm_r = sqrt(r_r);

norm_r0 = norm_r;

共軛梯度法初始時刻的共軛方向 并在代碼中以 mdelta 表示 (minus delta)。mdelta = r;

d_Pd = r_r;

e_Pd = 0;

共軛梯度法優化的目標函數。model_fun = @(eta, Heta) eta'*grad + .5*(eta'*Heta);

model_value = 0;

當迭代以達到最大迭代步數停止時，記 stop_tCG=5。stop_tCG = 5;

迭代主循環

最大迭代步數為 maxitfor j = 1 : opts.maxit

表示梯度， 表示海瑟矩陣。此處計算 。注意到這一步計算的耗時相對較長。Hmdelta = hess(mdelta);

計算曲率 和共軛梯度法的步長 。

在當前的搜索方向和步長下，我構造新的共軛方向 。計算 ，以便進行截斷共軛梯度法的檢查。為了簡化計算這里利用了 。d_Hd = mdelta'*Hmdelta;

alpha = r_r/d_Hd;

e_Pe_new = e_Pe + 2.0*alpha*e_Pd + alpha*alpha*d_Pd;

截斷共軛梯度法需要對曲率進行檢查，如果曲率滿足 ， 說明海瑟矩陣不正定(再求解下去得到的方向不一定是下降方向)，則停止共軛梯度算法。當 時， 迭代點超出可行域邊界，則通過構造得到一個位于邊界上的近似解并停止算法。

當上述兩條件中的任何一個成立，計算 使得 。 以 為最終的迭代結果。更新 。if d_Hd <= 0 || e_Pe_new >= Delta^2

tau = (-e_Pd + sqrt(e_Pd*e_Pd + d_Pd*(Delta^2-e_Pe))) / d_Pd;

eta = eta - tau*mdelta;

Heta = Heta -tau*Hmdelta;

記錄共軛梯度法的退出條件，當以非正曲率退出時記為 1，以超出邊界退出時記為 2。退出迭代。if d_Hd <= 0

stop_tCG = 1;

else

stop_tCG = 2;

end

break;

end

如果一步迭代更新的 沒有達到截斷共軛梯度法的兩個新增的停止條件，則更新 。e_Pe = e_Pe_new;

new_eta = eta-alpha*mdelta;

new_Heta = Heta -alpha*Hmdelta;

計算 處的目標函數值。如果一步更新后的目標函數值沒有下降，退出迭代，拒絕該步更新。 stop_tCG==6 表示目標函數值非降而終止。new_model_value = model_fun(new_eta, new_Heta);

if new_model_value >= model_value

stop_tCG = 6;

break;

end

如果沒有達到上述兩種情況，接受更新的 ，利用新的 更新變量。eta = new_eta;

Heta = new_Heta;

model_value = new_model_value;

更新殘差 以及其范數。特別的，記錄上一步的 為 r_rold。r = r -alpha*Hmdelta;

r_rold = r_r;

r_r = r'*r;

norm_r = sqrt(r_r);

標準的共軛梯度法的收斂條件：當達到最小迭代次數，且 時，認為算法收斂，停止迭代。

滿足上述收斂條件時，如果 ，則說明條件 為更嚴格的條件，說明此時外層牛頓法或者信賴域方法處于線性收斂的階段。 反之，條件 更嚴格，說明此時 已經較小，此時對應外層牛頓法或者信賴域方法的超線性收斂階段。if j >= opts.minit && norm_r <= norm_r0*min(norm_r0^theta, kappa)

if kappa < norm_r0^theta

stop_tCG = 3;

else

stop_tCG = 4;

end

break;

end

計算新的搜索方向： , 。beta = r_r/r_rold;

mdelta = r + beta*mdelta;

更新 ，這是由于 且 為 的線性組合，則由 可知 。

以及 ，注意到 ，有 。e_Pd = beta*(e_Pd + alpha*d_Pd);

d_Pd = r_r + beta*beta*d_Pd;

end

退出循環，記錄退出信息。out.iter = j;

out.stop_tCG = stop_tCG;end

參考頁面

此頁面實現了截斷共軛梯度算法，該函數用于 非精確牛頓法 中的牛頓方程的近似求解以及 信賴域方法 中的信賴域子問題的求解。

此頁面的源代碼請見： tCG.m。

版權聲明

此頁面為《最優化：建模、算法與理論》、《最優化計算方法》配套代碼。 代碼作者：文再文、劉浩洋、戶將，代碼整理與頁面制作：楊昊桐。

著作權所有 (C) 2020 文再文、劉浩洋、戶將

總結

以上是生活随笔為你收集整理的matlab cg steihaug,截断共轭梯度法的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。