Scipy优化Hawkes过程参数时直接返回初始值问题求助

我正在针对Hawkes过程开展随机模拟研究：

• 拥有n维真实向量D_True，由5个参数X={x1,x2,x3,x4,x5}驱动采样过程，得到n维向量D_Sim
• 得分计算方式：将D_True与D_Sim分别排序后，计算逐元素绝对差的均值
• 参数X需满足非线性约束，确保采样得到的alpha矩阵稳定，因此采用Scipy的SLSQP方法（L-BFGS-B无法处理约束）

我的代码如下：

import numpy as np

from tick.hawkes import SimuHawkesExpKernels
from scipy.sparse import rand
from scipy.spatial.distance import cityblock

n_nodes = 10  #dimension of vector D
D_true  = sorted[1797, 495, 7167, 4570, 8482, 5547, 524, 551, 525, 3490]) 

def distance(D_true, D_sampled):
    return 1/len(D_true) * cityblock(D_true, D_hat)

def simulate_alpha(n_nodes, k_r, theta_r, k_u, theta_u, gamma):
    
    rng = np.random.default_rng(seed = 332)

    u = rng.gamma(k_u, scale = 1/theta_u, size=(n_nodes, n_nodes))
    r =  rng.gamma(k_r, scale = 1/theta_r, size=(n_nodes, n_nodes))

    beta = gamma + u
    alpha = beta + r
    return alpha

def simulate_D(alpha, n_nodes = 10, end_time = 1000):
    decays = 5 * np.ones((n_nodes, n_nodes))
    baseline = 0.5 * np.ones(n_nodes)
    hawkes = SimuHawkesExpKernels(adjacency=alpha, decays=decays,
                                  baseline=baseline, verbose=False, end_time = end_time, seed=2398)
    dt = 0.1
    hawkes.track_intensity(dt)
    hawkes.simulate()
    
    simulated_D = = sorted([len(hawkes.timestamps[i]) for i in range(hawkes.n_nodes)])
    return simulated_D

score = lambda x: distance(D_true, 
                           simulate_D(simulate_alpha(n_nodes, x[0], x[1], x[2], 
                           x[3], x[4]))

stability_eq = lambda x: abs(eigvals(simulate_alpha(n_nodes=n_nodes,
                                                    k_r=x[0],
                                                    theta_r=x[1],
                                                    k_u=x[2],
                                                    theta_u=x[3],
                                                    gamma=x[4]))).max()

stability_constr = NonlinearConstraint(stability_eq, 0, 1)

x0 = [0.01, 0.06 * n_nodes, 0.01, 0.06 * n_nodes, 1 / (4*n_nodes)]

res = minimize(fun = score,
               x0=x0,
               method='SLSQP',
               constraints=stability_constraint,
               tol = 10)

问题现象 #

运行代码时，优化器仅完成1次迭代就返回初始值x0，此时得分仍为2250.5，并未收敛。scipy.minimize的输出如下：

fun: 2250.5
     jac: array([0., 0., 0., 0., 0.])
 message: 'Optimization terminated successfully'
    nfev: 6
     nit: 1
    njev: 1
  status: 0
 success: True
       x: array([0.01 , 0.6  , 0.01 , 0.6  , 0.025])

我尝试过其他求解器，但因优化问题非光滑，仅能使用SLSQP和COBYLA；同时尝试过不同初始值x0，结果仍返回初始值。

问题根源分析 #

1. 梯度异常：输出中jac: array([0., 0., 0., 0., 0.])是核心问题——SLSQP判定当前点梯度为0，直接终止。这是因为目标函数和约束是随机模拟驱动的非光滑、非连续函数：
   • simulate_alpha的随机采样导致函数输出对参数X的导数不存在，数值梯度噪声极大
   • 排序操作完全非光滑，进一步破坏梯度连续性
2. 代码语法错误：多处语法问题导致计算逻辑异常，间接影响优化：
   • D_true = sorted[1797, ...]应为D_true = sorted([1797, ...])
   • distance函数中未定义的D_hat应改为D_sampled
   • simulated_D = = sorted(...)多了一个等号
   • score的lambda表达式缺少闭合括号
   • 未导入scipy.optimize的minimize/NonlinearConstraint、numpy.linalg的eigvals

解决建议 #

1. 先修复代码语法错误

修正所有语法问题，确保函数能正确计算得分和约束值，这是优化的前提。

2. 处理随机模拟的梯度问题

由于目标函数是随机的，数值梯度不可靠，建议：

• 多次模拟求平均：对每组参数X，多次运行simulate_D取平均得分，降低噪声，让目标函数更平滑
• 切换到无梯度优化器：COBYLA更适合非光滑、随机目标，需调整参数：
  • 增大maxiter（比如设为1000）
  • 减小tol（当前10过大，建议设为1e-3）
  • 给参数设置正边界约束（gamma分布参数必须为正）

3. 改进约束计算稳定性

当前stability_eq每次重新采样alpha矩阵，约束值随机波动，优化器无法稳定判断：

• 改为计算alpha矩阵的期望谱半径：gamma分布期望为k/theta，因此E[alpha] = gamma + k_u/theta_u + k_r/theta_r，直接计算期望矩阵的谱半径作为约束，避免随机波动。

4. 调整优化参数设置

• 对SLSQP：若坚持使用，需手动提供近似梯度（或设置更小的eps步长），但随机函数梯度估计难度大
• 对COBYLA：添加正边界约束，示例配置：

  bounds = [(1e-5, None)] * 5
  res = minimize(fun=score, x0=x0, method='COBYLA',
                 constraints=stability_constr, tol=1e-3,
                 maxiter=1000, bounds=bounds)
  

修正后的核心代码片段

import numpy as np
from numpy.linalg import eigvals
from tick.hawkes import SimuHawkesExpKernels
from scipy.spatial.distance import cityblock
from scipy.optimize import minimize, NonlinearConstraint

n_nodes = 10  # dimension of vector D
D_true = sorted([1797, 495, 7167, 4570, 8482, 5547, 524, 551, 525, 3490]) 

def distance(D_true, D_sampled):
    return 1/len(D_true) * cityblock(D_true, D_sampled)

def simulate_alpha(n_nodes, k_r, theta_r, k_u, theta_u, gamma):
    rng = np.random.default_rng(seed=332)
    u = rng.gamma(k_u, scale=1/theta_u, size=(n_nodes, n_nodes))
    r = rng.gamma(k_r, scale=1/theta_r, size=(n_nodes, n_nodes))
    beta = gamma + u
    alpha = beta + r
    return alpha

def simulate_D(alpha, n_nodes=10, end_time=1000):
    decays = 5 * np.ones((n_nodes, n_nodes))
    baseline = 0.5 * np.ones(n_nodes)
    hawkes = SimuHawkesExpKernels(adjacency=alpha, decays=decays,
                                  baseline=baseline, verbose=False, end_time=end_time, seed=2398)
    hawkes.simulate()
    simulated_D = sorted([len(hawkes.timestamps[i]) for i in range(hawkes.n_nodes)])
    return simulated_D

# 多次模拟求平均得分，降低噪声
def score(x, n_sim=5):
    total_score = 0
    for _ in range(n_sim):
        alpha = simulate_alpha(n_nodes, x[0], x[1], x[2], x[3], x[4])
        D_sim = simulate_D(alpha)
        total_score += distance(D_true, D_sim)
    return total_score / n_sim

# 使用期望alpha矩阵计算谱半径，避免随机波动
def stability_eq(x):
    exp_u = x[2] / x[3]
    exp_r = x[0] / x[1]
    exp_alpha = x[4] + exp_u + exp_r
    exp_alpha_matrix = np.full((n_nodes, n_nodes), exp_alpha)
    return abs(eigvals(exp_alpha_matrix)).max()

stability_constr = NonlinearConstraint(stability_eq, 0, 1)
x0 = [0.01, 0.06 * n_nodes, 0.01, 0.06 * n_nodes, 1 / (4*n_nodes)]
bounds = [(1e-5, None)] * 5  # 所有参数必须为正

# 使用COBYLA尝试优化
res = minimize(fun=score, x0=x0, method='COBYLA',
               constraints=stability_constr, tol=1e-3,
               maxiter=1000, bounds=bounds)
print(res)

内容的提问来源于stack exchange，提问作者kohlstein