使用SLSQP求解最大夏普比率组合时出现分母NaN问题求助

问题场景 #

使用Scipy的SLSQP算法求解最大夏普比率投资组合时，目标函数的分母始终输出NaN（自定义回调函数可验证），但将初始权重x0代入分母计算能得到有效数值，疑惑SLSQP运行初期为何出现该问题。

用户原代码如下：

from scipy.optimize import minimize
import numpy as np

n = 250
x0 = np.ones(n) / n
pred = np.random.normal(loc=0.07, scale=0.2, size=(n))
c = ((np.random.random(size = (n, n)))-0.5) *0.1

def objective(weights):
  if np.sum(weights) == 0:
    return 10
  else:
    return -np.dot(weights, pred)/(np.sqrt(np.dot(weights, np.dot(c, weights))))

def custom_callback(xk):
    print(np.dot(xk, pred) , (np.sqrt(np.dot(xk, np.dot(c, xk)))))

constraints = (
  {'type': 'ineq', 'fun': lambda weights: np.abs(np.sum(weights)) - 0.03},  
  {'type': 'ineq', 'fun': lambda weights: 1.03 - np.sum(np.abs(weights))},  
  {'type': 'ineq', 'fun': lambda weights: np.sum(np.abs(weights))-0.97},
  )

bounds = tuple((-0.02, 0.02) for asset in range(len(x0)))
result = minimize(objective, x0, method='SLSQP', bounds=bounds, constraints=constraints, options={'disp': 2, 'maxiter' : 100}, callback=custom_callback)

验证初始权重分母的代码：

np.sqrt(np.dot(x0, np.dot(c, x0)))

问题根源 #

1. 协方差矩阵非半正定：投资组合的协方差矩阵必须是半正定的（二次型结果≥0），但用户生成c的方式((np.random.random(size=(n,n)))-0.5)*0.1会生成随机非对称矩阵，且大概率不满足半正定条件。当SLSQP迭代过程中尝试某些权重时，二次型np.dot(weights, np.dot(c, weights))可能得到负数，开根号后就会产生NaN。
2. 初始点巧合：初始权重x0恰好让二次型结果为正，所以计算正常，但迭代中的权重组合没有这个保证。

解决方案 #

1. 生成合法的半正定协方差矩阵 #

将协方差矩阵c改为对称半正定矩阵，比如通过随机矩阵的转置乘积生成，再添加小正则项避免数值奇异：

# 生成对称半正定协方差矩阵
a = np.random.randn(n, n) * 0.1
c = np.dot(a.T, a) + 1e-6 * np.eye(n)  # 添加小对角项保证正定

2. 目标函数添加数值保护 #

在计算分母前先检查二次型结果，若小于等于阈值则返回惩罚值，避免NaN：

def objective(weights):
    weight_sum = np.sum(weights)
    if weight_sum == 0:
        return 10
    var = np.dot(weights, np.dot(c, weights))
    # 用小阈值处理数值误差，避免开根号出现NaN
    if var <= 1e-10:
        return 10  # 返回大的惩罚值，引导优化器避开该点
    return -np.dot(weights, pred) / np.sqrt(var)

3. 优化回调函数 #

在回调函数中也添加检查，避免输出NaN：

def custom_callback(xk):
    mean = np.dot(xk, pred)
    var = np.dot(xk, np.dot(c, xk))
    std = np.sqrt(var) if var > 1e-10 else np.nan
    print(mean, std)

修改后完整代码 #

from scipy.optimize import minimize
import numpy as np

n = 250
x0 = np.ones(n) / n
pred = np.random.normal(loc=0.07, scale=0.2, size=(n))
# 生成合法的对称半正定协方差矩阵
a = np.random.randn(n, n) * 0.1
c = np.dot(a.T, a) + 1e-6 * np.eye(n)

def objective(weights):
    weight_sum = np.sum(weights)
    if weight_sum == 0:
        return 10
    var = np.dot(weights, np.dot(c, weights))
    if var <= 1e-10:
        return 10
    return -np.dot(weights, pred) / np.sqrt(var)

def custom_callback(xk):
    mean = np.dot(xk, pred)
    var = np.dot(xk, np.dot(c, xk))
    std = np.sqrt(var) if var > 1e-10 else np.nan
    print(mean, std)

constraints = (
  {'type': 'ineq', 'fun': lambda weights: np.abs(np.sum(weights)) - 0.03},  
  {'type': 'ineq', 'fun': lambda weights: 1.03 - np.sum(np.abs(weights))},  
  {'type': 'ineq', 'fun': lambda weights: np.sum(np.abs(weights))-0.97},
  )

bounds = tuple((-0.02, 0.02) for asset in range(len(x0)))
result = minimize(objective, x0, method='SLSQP', bounds=bounds, constraints=constraints, options={'disp': 2, 'maxiter' : 100}, callback=custom_callback)

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