如何在scipy差分进化优化中验证方差-协方差矩阵且不影响收敛？

我正在开展一项以方差-协方差矩阵为待估计参数的优化任务，使用scipy库的differential_evolution工具实现优化。优化前100步中，回调显示f(x)=inf，经验证这些步长对应的参数均不满足约束条件。但这种情况导致了收敛问题：实际任务中已运行458步，convergence仍为0.0。现咨询如何在不影响优化收敛的前提下，验证方差-协方差矩阵的有效性？

原实现代码：

import numpy as np
from scipy.optimize import differential_evolution, NonlinearConstraint

def likelihood(params):
    p = params[0] # this parameter is not related to var-covar matrix
    dev = np.diag(params[8:15])
    X = np.eye(dev.shape[0])
    X[0, 1:7] = params[15:21]
    X[1, 2:7] = params[21:26]
    X[2, 3:7] = params[26:30]
    X[3, 4:7] = params[30:33]
    X[4, 5:7] = params[33:35]
    X[5, 6:7] = params[35:36]
    X = X + X.T - np.diag(np.diag(X))
    cov = dev.dot(X).dot(dev)
    L = np.linalg.cholesky(cov)
    
    # only for reproduce.
    arg1 = -0.2 * np.sqrt(0.5 * (params[0] ** 2 + params[1] ** 2))
    arg2 = 0.5 * (np.cos(2. * np.pi * params[0]) + np.cos(2. * np.pi * params[1]))
    return -20. * np.exp(arg1) - np.exp(arg2) + 20. + np.e + params.sum()
    

def constrain(params):
    dev = np.diag(params[8:15])
    X = np.eye(dev.shape[0])
    X[0, 1:7] = params[15:21]
    X[1, 2:7] = params[21:26]
    X[2, 3:7] = params[26:30]
    X[3, 4:7] = params[30:33]
    X[4, 5:7] = params[33:35]
    X[5, 6:7] = params[35:36]
    X = X + X.T - np.diag(np.diag(X))
    cov = dev.dot(X).dot(dev)
    try:
        np.linalg.cholesky(cov)
        return 0
    except np.linalg.LinAlgError:
        return 1

def print_de(x, convergence):
    print(x,"convergence: ",convergence)

bounds = (((0.0, 1.0),) * 8 + ((0.0, 1.0),) + ((0.0, 0.5),) * 6 + ((-1.0,1.0),)*21)
differential_evolution(likelihood, bounds=bounds, disp=True,
                       callback=print_de,
                       constraints=NonlinearConstraint(constrain, 0, 0)
                      )

1. 重构协方差矩阵生成逻辑，消除冗余代码 #

先把重复的协方差矩阵生成逻辑提取为独立函数，减少计算开销，也方便后续修改：

def build_covariance(params):
    dev = np.diag(params[8:15])
    X = np.eye(dev.shape[0])
    X[0, 1:7] = params[15:21]
    X[1, 2:7] = params[21:26]
    X[2, 3:7] = params[26:30]
    X[3, 4:7] = params[30:33]
    X[4, 5:7] = params[33:35]
    X[5, 6:7] = params[35:36]
    X = X + X.T - np.diag(np.diag(X))
    return dev.dot(X).dot(dev)

2. 替换事后约束检查为天然正定的参数化方式（最优方案） #

原方案通过非线性约束强制协方差矩阵正定，效率低且容易产生大量无效解。更高效的方式是直接参数化正定矩阵：用下三角矩阵L（对角线元素为正）构造cov = L @ L.T，这样生成的矩阵天然正定，无需额外约束检查。

修改参数映射逻辑：

• 原参数中对应协方差的部分（params[8:]）替换为下三角矩阵L的元素（对角线设为正区间，非对角线设为任意实数区间）
• 重构协方差矩阵生成函数为基于下三角矩阵的方式

示例修改后的代码：

def build_covariance_from_lower_tri(params_lower_tri):
    # params_lower_tri包含下三角矩阵的元素，长度为7*(7+1)/2=28
    L = np.zeros((7,7))
    idx = 0
    for i in range(7):
        # 对角线元素设为正区间，保证矩阵正定
        L[i,i] = params_lower_tri[idx]
        idx +=1
        for j in range(i):
            L[i,j] = params_lower_tri[idx]
            idx +=1
    return L @ L.T

# 重新定义边界：前8个参数保持原边界，后28个参数对应下三角矩阵
# 对角线元素（7个）设为(0.01, 1.0)避免接近0，非对角线21个设为(-1.0,1.0)
bounds = (((0.0, 1.0),)*8) + (((0.01, 1.0),)*7) + (((-1.0,1.0),)*21)

def likelihood(params):
    p = params[0]
    # 提取下三角矩阵参数
    cov_params = params[8:]
    cov = build_covariance_from_lower_tri(cov_params)
    # 原目标函数逻辑保留
    arg1 = -0.2 * np.sqrt(0.5 * (params[0] ** 2 + params[1] ** 2))
    arg2 = 0.5 * (np.cos(2. * np.pi * params[0]) + np.cos(2. * np.pi * params[1]))
    return -20. * np.exp(arg1) - np.exp(arg2) + 20. + np.e + params.sum()

这种方式彻底消除了无效解，种群所有个体的参数都能生成合法的协方差矩阵，大幅提升收敛效率。

3. 若必须保留原参数化，优化约束与目标函数处理 #

如果无法修改参数化方式，可通过以下方式优化：

• 在目标函数中，当协方差矩阵不正定时，返回一个极大的惩罚值（而非让Cholesky报错导致inf），避免种群陷入无效解区域
• 简化约束函数，利用矩阵正定的等价条件（所有特征值为正）或更高效的检查方式

修改后的目标函数与约束：

def likelihood(params):
    p = params[0]
    cov = build_covariance(params)
    # 快速检查正定：计算最小特征值，若小于等于0则返回惩罚值
    min_eig = np.linalg.eigvalsh(cov).min()
    if min_eig <= 1e-8:
        # 返回一个极大值作为惩罚，避免inf
        return 1e10 + params.sum()
    # 原目标函数逻辑
    arg1 = -0.2 * np.sqrt(0.5 * (params[0] ** 2 + params[1] ** 2))
    arg2 = 0.5 * (np.cos(2. * np.pi * params[0]) + np.cos(2. * np.pi * params[1]))
    return -20. * np.exp(arg1) - np.exp(arg2) + 20. + np.e + params.sum()

def constrain(params):
    cov = build_covariance(params)
    # 用特征值检查替代Cholesky，更稳定
    min_eig = np.linalg.eigvalsh(cov).min()
    # 返回最小特征值，约束其>=0（留微小epsilon避免数值误差）
    return min_eig - 1e-8

同时修改约束定义：

constraints = NonlinearConstraint(constrain, 0, np.inf)

这种方式减少了inf的出现，让优化器能更平滑地引导种群向可行区域移动。

4. 调整differential_evolution的参数加速收敛 #

• 增大popsize（比如设置为popsize=20），增加种群多样性，更快探索可行区域
• 调整mutation和recombination参数，平衡探索与 exploitation
• 启用polish=True，在优化后期用局部搜索精化解

示例：

differential_evolution(likelihood, bounds=bounds, disp=True,
                       callback=print_de,
                       constraints=constraints,
                       popsize=20,
                       mutation=(0.5, 1.0),
                       recombination=0.7,
                       polish=True
                      )

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