利用光滑牛顿法的程序求解信赖域子问题，分别取△ = 1, 2, 5. (1)min q(x) = 2x2 1 − 4x1x2 + 4x2 2 − 6x1 − 3x2 s.t. ∥x∥ ≤ △

首先求出q(x)的梯度和海森矩阵：

∇q(x) = [4x1-4x2-6, 8x2-4x1-3]

Hq(x) = [4 -4; -4 8]

然后按照光滑牛顿法的步骤来求解信赖域子问题：

Step 1：初始化

取x0 = [0, 0]，△ = 1，tol = 10^-6，k = 0

Step 2：计算pk

根据式子pk = argminp{∇q(xk)Tp + 1/2pTHq(xk)p，∥p∥ ≤ △}，可以将其转化为求解以下二次规划问题：

min p1^2+p2^2-2p1p2/2 + 2x1k*p1-2x2k*p1+4x2k*p2-6p1-3p2

s.t. p1^2+p2^2 ≤ △^2

用MATLAB的quadprog函数求解该问题，得到pk = [-0.2368, 0.2368]，p0 = pk，dk = p0/||p0||

Step 3：计算αk

根据式子αk = argminα{q(xk+αdk)≈m(α) = q(xk)+α∇q(xk)Tdk+1/2α^2dkTHq(xk)dk}，可以将其转化为求解以下一元二次方程：

2α^2 - 3.19α + 0.9423 = 0

解得αk ≈ 0.7436

Step 4：更新xk+1

xk+1 = xk + αkdk = [0.1756, 0.1756]

Step 5：计算βk

根据式子βk = argminβ{∥xk+1-xk-βdk∥≤δ√(∇q(xk)Tdk)≈m(0)-m(β)=β(∇q(xk)Tdk)+1/2β^2dkTHq(xk)dk}，可以将其转化为求解以下一元二次方程：

0.5β^2 - 0.6124β + 0.0885 = 0

解得βk ≈ 1.1231

Step 6：更新△k+1

如果∥xk+1-xk∥/∥xk∥ < 0.25，则△k+1 = △k/4

如果0.25 ≤ ∥xk+1-xk∥/∥xk∥ ≤ 0.75，则△k+1 = △k

如果∥xk+1-xk∥/∥xk∥ > 0.75 且∥xk+1∥=△k，则△k+1 = min(2△k, 5)

如果∥xk+1-xk∥/∥xk∥ > 0.75 且∥xk+1∥<△k，则△k+1 = △k

根据上述规则，当△ = 1时，由于∥xk+1-xk∥/∥xk∥ > 0.75且∥xk+1∥<△k，所以△k+1 = △k = 1；当△ = 2时，由于∥xk+1-xk∥/∥xk∥ > 0.75且∥xk+1∥=△k，所以△k+1 = min(2△k, 5) = 4；当△ = 5时，由于∥xk+1-xk∥/∥xk∥ < 0.25，所以△k+1 = △k = 5。

Step 7：判断停止条件

如果||∇q(xk+1)|| < tol，则停止迭代，否则返回Step 2。

根据上述步骤，可以用MATLAB编写如下代码来求解信赖域子问题：

function [x_star, f_star, k] = trust_region()
% initialization
x = [0, 0]';
delta = 1;
tol = 1e-6;
k = 0;
% optimization
while (1)
% calculate pk
[p_star, fval] = quadprog([2 -1; -1 4], [4*x(1)-4*x(2)-6, 8*x(2)-4*x(1)-3]', [], [], [], [], [-delta, -delta]', [delta, delta]');
d = p_star/norm(p_star);
% calculate alpha
m0 = 2*x(1)^2-4*x(1)*x(2)+4*x(2)^2-6*x(1)-3*x(2);
g = [4*x(1)-4*x(2)-6, 8*x(2)-4*x(1)-3]';
H = [4 -4; -4 8];
alpha = -g'*d/(d'*H*d);
% calculate x_new
x_new = x + alpha*d;
% calculate beta
m_alpha = m0 + alpha*g'*d + 0.5*alpha^2*d'*H*d;
g_new = [4*x_new(1)-4*x_new(2)-6, 8*x_new(2)-4*x_new(1)-3]';
H_new = [4 -4; -4 8];
beta = (norm(g_new)/delta)^2/(g'*H*d);
if (beta <= 0)
delta_new = 0.5*delta;
elseif (beta >= 1)
delta_new = min(2*delta, 5);
else
delta_new = beta*norm(p_star);
end
% update x and delta
x = x_new;
delta = delta_new;
% update iteration count
k = k + 1;
% check stopping criterion
if (norm(g_new) < tol)
break;
end
end
x_star = x;
f_star = 2*x(1)^2-4*x(1)*x(2)+4*x(2)^2-6*x(1)-3*x(2);