Numpy矩阵乘法语法及有效前沿双曲线代码问题排查

我想用Numpy把给定公式转成Python代码，其中mu是3×1向量，sigma是3×3协方差矩阵。最初写的代码输出了3×3×3的结果，调整结构后，画出的曲线还是和预期的红色曲线不符，麻烦排查问题，包括向量转置、矩阵乘法这类可能的错误。

初始代码 #

import numpy as np

u = np.array([0.0493, 0.0770, 0.0886])
cov_matrix = ([[0.0017, -0.0017, -0.0021],
               [-0.0017, 0.0396, 0.0309],
               [-0.0021, 0.0309, 0.0392]])

def hyperbola(u):
    k = u * np.linalg.inv(cov_matrix)                       
    l = u * np.linalg.inv(cov_matrix) * u[:, np.newaxis]    
    m = np.linalg.inv(cov_matrix)                           

    g = (l * np.linalg.inv(cov_matrix) - k * np.linalg.inv(cov_matrix) * u[:, np.newaxis]) / (l*m - k**2)
    h = (m * np.linalg.inv(cov_matrix) * u[:, np.newaxis] - k * np.linalg.inv(cov_matrix)) / (l*m - k**2)
    
    a = h * cov_matrix * h[:, np.newaxis]
    b = 2 * g * cov_matrix * h[:, np.newaxis]
    c = g * cov_matrix * g[:, np.newaxis]

    return np.sqrt(a * u**2 + b*u + c)


efficient_frontier = hyperbola(u)

调整后代码 #

import numpy as np
import matplotlib.pyplot as plt

u = np.array([0.0493, 0.0770, 0.0886])
cov_matrix = np.array([[0.0017, -0.0017, -0.0021],
                       [-0.0017, 0.0396, 0.0309],
                       [-0.0021, 0.0309, 0.0392]])

def hyperbola(u):
    k = u @ np.linalg.inv(cov_matrix) @ np.ones_like(u)                  
    l = u @ np.linalg.inv(cov_matrix) @ u
    m = np.ones_like(u) @ np.linalg.inv(cov_matrix) @ np.ones_like(u)

    g = (l * np.linalg.inv(cov_matrix) @ np.ones_like(u) - k * np.linalg.inv(cov_matrix) @ u) / (l*m - k**2)
    h = (m * np.linalg.inv(cov_matrix) @ u - k * np.linalg.inv(cov_matrix) @ np.ones_like(u)) / (l*m - k**2)
    
    a = h @ cov_matrix @ h
    b = 2 * g @ cov_matrix @ h
    c = g @ cov_matrix @ g

    return np.sqrt(a * u**2 + b*u + c)

efficient_frontier = hyperbola(u)

plt.plot(u, efficient_frontier)
plt.xlabel('Volatility')
plt.ylabel('Mean Return')
plt.grid(True)

问题分析与修正 #

核心错误点 #

1. 矩阵乘法误用：初始代码用*做元素级乘法，而非@矩阵乘法，直接导致维度爆炸（3×3×3结果）。调整后的代码虽改用@，但g、h的计算逻辑仍有误——标量与矩阵相乘后再和向量运算，维度匹配但逻辑不符合有效前沿公式。
2. 有效前沿逻辑误解：原代码直接用输入的u（资产期望收益向量）计算输出，这是单个资产的收益-波动率对应，而非遍历连续目标收益得到的有效前沿曲线。有效前沿需要生成一系列目标收益值，计算对应最小波动率。
3. 向量维度处理：g、h应为行向量，后续与协方差矩阵的乘法需严格遵循矩阵维度规则，原代码中部分乘法会导致维度不匹配或结果错误。

修正后代码 #

import numpy as np
import matplotlib.pyplot as plt

# 输入参数
mu = np.array([0.0493, 0.0770, 0.0886])
cov_matrix = np.array([[0.0017, -0.0017, -0.0021],
                       [-0.0017, 0.0396, 0.0309],
                       [-0.0021, 0.0309, 0.0392]])

# 预计算逆矩阵与常量
inv_cov = np.linalg.inv(cov_matrix)
ones = np.ones_like(mu)

k = mu.T @ inv_cov @ ones
l = mu.T @ inv_cov @ mu
m = ones.T @ inv_cov @ ones

# 计算g、h行向量
g = (l * (ones.T @ inv_cov) - k * (mu.T @ inv_cov)) / (l * m - k**2)
h = (m * (mu.T @ inv_cov) - k * (ones.T @ inv_cov)) / (l * m - k**2)

# 计算a、b、c标量
a = h @ cov_matrix @ h.T
b = 2 * g @ cov_matrix @ h.T
c = g @ cov_matrix @ g.T

# 有效前沿函数：输入目标收益数组，输出对应最小波动率
def efficient_frontier(mu_target):
    sigma_sq = a * mu_target**2 + b * mu_target + c
    return np.sqrt(np.maximum(sigma_sq, 0))  # 确保波动率非负

# 生成连续目标收益范围
mu_range = np.linspace(mu.min() - 0.02, mu.max() + 0.02, 100)
sigma_range = efficient_frontier(mu_range)

# 绘图：有效前沿+单个资产点
plt.plot(sigma_range, mu_range, label='有效前沿')
asset_vols = np.sqrt(np.diag(cov_matrix))
plt.scatter(asset_vols, mu, color='red', label='单个资产')
plt.xlabel('波动率')
plt.ylabel('期望收益')
plt.grid(True)
plt.legend()
plt.show()

修正说明 #

• 预计算逆矩阵与所有常量，避免重复计算提升效率
• 调整g、h的计算逻辑，确保符合均值-方差有效前沿的数学公式
• 重新定义有效前沿函数，接收连续目标收益数组，输出对应波动率
• 添加单个资产的红色散点，与有效前沿曲线做对比，验证结果正确性

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