如何利用信赖域方法从鞍点引导优化至指定极小值点A？

现有基于Jorge Nocedal & Stephen Wright《数值优化》（第二版）实现的信赖域算法，从鞍点启动时会收敛至极小值点C，需调整算法使其从鞍点出发后收敛至预先指定的极小值点A。原核心代码如下：

''' ----------------------------------------------------------------------
Trust-Region Methods Jorge Nocedal & Stephen Wright. 
Numerical Optimization, second edition.
--------------------------------------------------------------------------'''
def trm(xk, alfa, deltak, maxiter=100, maxdelta=100.0, etha=0.25, tol=1e-1):
    k = 0
    xk0 = []
    xk1 = []
    fks = []
    ks = []
    while True:
        fk = func(xk)
        gk = grad(xk)
        Bk = hess(xk)
        pk, deltak = subproblem(gk, Bk, deltak)
        
        rhok = (fk - func(xk + pk)) / -(np.dot(gk.T, pk) + 0.5 * np.dot(np.dot(pk.T, Bk), pk))
        if rhok < 0.25:
            deltak = 0.25 * deltak
        else:
            if rhok < 0.75 and la.norm(pk) == deltak:
                deltak = min(2* deltak, maxdelta)
            else:
                deltak = deltak
        if rhok > etha:
            xk = xk + pk * alfa
        else:
            xk = xk
        
        xk0.append(float(xk[0]))
        xk1.append(float(xk[1]))

        fks.append(float(fk))
        ks.append(int(k))

        k += 1
        if la.norm(gk) < tol or k > maxiter:
            print("End point:", xk)
            print("Energy:", fk)
            plot_energy(ks, fks)
            plot_countour(xk0, xk1)
            plt.show()
            break
    
    return xk0, xk1


def subproblem(g, B, delta):
    b, Q = la.eigh(B)
    f = np.dot(Q.T, g)
    lam = -min(b)
    if abs(f[0]) > CLOSE_TO_ZERO:
        phi = lambda x: (1 / delta) - (1 / np.sqrt(sum((f / (b + x))**2)))
        sol = root(phi, lam)
        lam = sol.x
        p = np.zeros((N_COORDS))
        for i in range(N_COORDS):
            p += -(f[i] / (b[i] + lam)) * Q[:, i]
        delta = la.norm(p)
    else:
        tau = 0
        for i in range(1, N_COORDS):
            tau += (f[i] / (b[i] - b[0]))**2
        tau = np.sqrt(abs(delta**2 - tau))
        p = np.zeros((N_COORDS))
        for i in range(1, N_COORDS):
            p += -(f[i] / (b[i] + lam)) * Q[:, i]
        delta = tau
    
    return p, delta

调整思路 #

鞍点处梯度为0，原算法的搜索方向由Hessian矩阵的负曲率方向主导，导致收敛到非目标极小值点。要定向收敛到A，需在信赖域的二次近似模型中加入目标点引导项，引导搜索方向同时兼顾局部二次拟合和向A靠近的趋势。

具体实现步骤 #

1. 新增目标引导参数与修改子问题模型 #

在子问题中引入指向目标点A的正则化项，修改二次近似模型：
原模型：$m_k(p) = f_k + g_k^T p + \frac{1}{2}p^T B_k p$
修改后：$m_k(p) = f_k + g_k^T p + \frac{1}{2}p^T B_k p + \lambda_g |x_k + p - A|^2$
其中$\lambda_g$是正权重参数，控制引导强度。

2. 重写子问题求解函数 #

更新subproblem函数，传入目标点A、当前点xk和引导权重，调整梯度与Hessian矩阵后求解：

def subproblem(g, B, delta, target_A, xk, lambda_g=1.0):
    # 加入引导项后的新梯度和Hessian
    new_g = g + 2 * lambda_g * (xk - target_A)
    new_B = B + 2 * lambda_g * np.eye(len(g))
    
    b, Q = la.eigh(new_B)
    f = np.dot(Q.T, new_g)
    lam = -min(b)
    CLOSE_TO_ZERO = 1e-8  # 需确保该常量已定义
    
    if abs(f[0]) > CLOSE_TO_ZERO:
        phi = lambda x: (1 / delta) - (1 / np.sqrt(sum((f / (b + x))**2)))
        sol = root(phi, lam)
        lam = sol.x
        p = np.zeros((len(g)))
        for i in range(len(g)):
            p += -(f[i] / (b[i] + lam)) * Q[:, i]
        delta = la.norm(p)
    else:
        tau = 0
        for i in range(1, len(g)):
            tau += (f[i] / (b[i] - b[0]))**2
        tau = np.sqrt(abs(delta**2 - tau))
        p = np.zeros((len(g)))
        for i in range(1, len(g)):
            p += -(f[i] / (b[i] + lam)) * Q[:, i]
        # 补充沿第一个特征向量的分量，满足信赖域约束
        p += tau * Q[:, 0]
        delta = la.norm(p)
    
    return p, delta

3. 更新主迭代函数 #

在trm函数中新增目标点参数，传递给子问题，并保留原信赖域的迭代逻辑：

def trm(xk, alfa, deltak, target_A, lambda_g=1.0, maxiter=100, maxdelta=100.0, etha=0.25, tol=1e-1):
    k = 0
    xk0 = []
    xk1 = []
    fks = []
    ks = []
    while True:
        fk = func(xk)
        gk = grad(xk)
        Bk = hess(xk)
        # 传入目标点与引导权重
        pk, deltak = subproblem(gk, Bk, deltak, target_A, xk, lambda_g)
        
        rhok = (fk - func(xk + pk)) / -(np.dot(gk.T, pk) + 0.5 * np.dot(np.dot(pk.T, Bk), pk))
        if rhok < 0.25:
            deltak = 0.25 * deltak
        else:
            if rhok < 0.75 and la.norm(pk) == deltak:
                deltak = min(2* deltak, maxdelta)
        if rhok > etha:
            xk = xk + pk * alfa
        
        xk0.append(float(xk[0]))
        xk1.append(float(xk[1]))
        fks.append(float(fk))
        ks.append(int(k))

        k += 1
        if la.norm(gk) < tol or k > maxiter:
            print("End point:", xk)
            print("Energy:", fk)
            plot_energy(ks, fks)
            plot_countour(xk0, xk1)
            plt.show()
            break
    
    return xk0, xk1

4. 参数调优建议 #

• lambda_g设置：初始值建议取1.0~10.0，若收敛方向偏离A则增大权重；若局部搜索能力下降则减小权重。
• 自适应调整：可根据当前点xk到A的距离动态调整lambda_g，距离越远权重越大，例如：

  lambda_g = 10.0 * la.norm(xk - target_A) / la.norm(target_A)
  

• 鞍点特殊处理：鞍点处梯度为0，引导项成为搜索方向的主导因素，确保第一步就朝向A移动。

调用示例 #

启动算法时传入指定的极小值点A：

# 假设A是二维点[1.0, 2.0]
target_A = np.array([1.0, 2.0])
# 从鞍点xk_saddle启动
xk_saddle = np.array([0.0, 0.0])
trm(xk_saddle, alfa=1.0, deltak=1.0, target_A=target_A, lambda_g=5.0)

内容的提问来源于stack exchange，提问作者Arturo Rentería