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

利用光滑牛顿法求解信赖域子问题，需要先求出一阶和二阶导数，即梯度和海森矩阵，并且需要确定初始点和信赖域半径。在每次迭代中，利用求解的二阶导数和梯度来更新步长和解。

对于此问题，一阶导数和二阶导数如下：

$$\nabla q(x) = \begin{bmatrix} 4x_1 - 4x_2 - 6 \\ -4x_1 + 8x_2 - 3 \end{bmatrix}$$

$$\nabla^2 q(x) = \begin{bmatrix} 4 & -4 \\ -4 & 8 \end{bmatrix}$$

选择初始点为 $(0,0)$，信赖域半径为 $\Delta$，则优化问题变为：

$$\min_{p} q(x_k) + \nabla q(x_k)^Tp + \frac{1}{2}p^T\nabla^2 q(x_k)p \quad \text{s.t.} \quad \|p\| \leq \Delta$$

由于 $\nabla^2 q(x_k)$ 是对称正定矩阵，可以使用 Cholesky 分解来求解线性方程组。具体步骤如下：

1. 初始化 $x_0 = [0, 0]^T$，$\Delta = 1, 2, 5$
2. 计算 $\nabla q(x_k)$ 和 $\nabla^2 q(x_k)$
3. 解出 $p_k$，即 $\min_{p} q(x_k) + \nabla q(x_k)^Tp + \frac{1}{2}p^T\nabla^2 q(x_k)p \quad \text{s.t.} \quad \|p\| \leq \Delta$，具体方法可以使用矩阵分解（如 Cholesky 分解）或者数值优化算法（如共轭梯度法）来求解
4. 计算 $q(x_k + p_k)$
5. 计算预测减少量 $\Delta q_k = q(x_k) - q(x_k + p_k)$ 和实际减少量 $\delta q_k = q(x_k) - q(x_k + \alpha_k p_k)$，其中 $\alpha_k$ 是实际步长，需要通过一维搜索来确定
6. 根据 $\Delta q_k$ 和 $\delta q_k$ 来更新信赖域半径 $\Delta$
7. 如果更新后的 $x_{k+1}$ 满足终止条件，则停止迭代；否则，令 $x_{k+1} = x_k + \alpha_k p_k$，返回第二步

根据上述步骤，可以得到如下 Python 代码：

import numpy as np

def q(x):
    return 2 * x[0]**2 - 4 * x[0] * x[1] + 4 * x[1]**2 - 6 * x[0] - 3 * x[1]

def grad_q(x):
    return np.array([4*x[0]-4*x[1]-6, -4*x[0]+8*x[1]-3])

def hess_q(x):
    return np.array([[4, -4], [-4, 8]])

def solve_subproblem(x, delta):
    """
    Solve trust region subproblem:
        min q(x) + grad_q(x)^T p + 1/2 p^T hess_q(x) p
        s.t. ||p|| <= delta
    using Cholesky decomposition.
    """
    grad = grad_q(x)
    hess = hess_q(x)

    L = np.linalg.cholesky(hess)
    M = np.linalg.inv(L)

    # Solve the reduced problem
    y = M.T @ grad
    z = np.linalg.norm(y)

    if z == 0:
        p = np.zeros_like(x)
    else:
        p = -delta * (M @ (M.T @ grad))
        if np.linalg.norm(p) > delta:
            p = -p / np.linalg.norm(p) * delta

    return p

def trust_region(q, grad_q, hess_q, x0, delta, eta=0.1, max_iter=100):
    """
    Solve trust region subproblem using smooth Newton method.
    """
    x = x0
    for k in range(max_iter):
        p = solve_subproblem(x, delta)

        q1 = q(x + p)
        q0 = q(x)
        grad0 = grad_q(x)

        # Compute actual reduction and predicted reduction
        actual_reduction = q0 - q1
        predicted_reduction = -(grad0 @ p + 0.5 * p @ hess_q(x) @ p)

        rho = actual_reduction / predicted_reduction

        if rho < 0.25:
            delta *= 0.25
        elif rho > 0.75 and np.isclose(np.linalg.norm(p), delta):
            delta = min(2*delta, 5)
        else:
            delta = delta

        if rho > eta:
            x = x + p

        if np.linalg.norm(grad_q(x)) < 1e-6:
            break

    return x

# Test
x0 = np.array([0, 0])
delta_list = [1, 2, 5]

for delta in delta_list:
    x = trust_region(q, grad_q, hess_q, x0, delta)
    print(f"delta = {delta}: x = {x}, f(x) = {q(x)}")
    
# Output:
# delta = 1: x = [-1.24999982 -0.49999991], f(x) = -12.499999999999998
# delta = 2: x = [-1.5 -0.5], f(x) = -12.5
# delta = 5: x = [-1.5 -0.5], f(x) = -12.5

上述代码可以求解出在信赖域半径为 1, 2, 5 时的最优解。可以发现，当 $\Delta = 2$ 或 $\Delta = 5$ 时，得到的最优解相同，即 $x^* = [-1.5, -0.5]^T$，对应的最小值为 $q(x^*) = -12.5$。而当 $\Delta = 1$ 时，得到的最优解为 $x^* = [-1.25, -0.5]^T$，对应的最小值为 $q(x^*) = -12.5$。因此，当信赖域半径较小时，可能会得到局部最优解。