1. 从经典难题到量子视角投资组合优化为何值得一试如果你在金融、量化或者运筹优化领域摸爬滚打过一阵子大概率听说过“投资组合优化”这个经典问题。简单来说就是给你一笔钱和一堆可供选择的资产比如股票、债券你需要在风险和收益之间找到一个最优的平衡点决定每样资产买多少。这听起来像是个简单的数学题但一旦资产数量上去约束条件复杂起来比如不允许卖空、有交易成本、行业配置限制它就会迅速变成一个计算上的“硬骨头”——一个典型的NP-hard组合优化问题。传统的解法比如二次规划对于中小规模问题还能应付但面对成百上千种资产和复杂的现实约束时要么求解时间长得无法接受要么干脆找不到最优解只能给个近似解。这也是为什么这个领域一直是学术界和工业界的前沿战场。而最近几年一个听起来有点科幻的技术——量子计算特别是其中的量子退火分支开始进入大家的视野。它提供了一种全新的思路将优化问题映射到物理系统的能量最低状态让量子比特通过量子隧穿效应“翻山越岭”去寻找那个全局最优解。你可能觉得这离我们太远需要昂贵的量子计算机。但有意思的是像D-Wave这样的公司已经提供了云端的量子退火服务而我们作为开发者完全可以用熟悉的Python工具链来构建和提交问题。这就是PyQUBO的用武之地。它是一个Python库能让你用声明式的方法定义优化问题目标函数和约束条件然后自动将其编译成QUBO二次无约束二进制优化形式这是量子退火机“听得懂”的语言。MathorCup2023将其作为赛题正是看中了它在解决复杂优化问题上的潜力和对前沿技术的实践价值。所以这篇内容不是飘在天上的理论探讨而是一次接地气的实战记录。我将带你一步步用PyQUBO把经典的投资组合优化问题“翻译”成QUBO模型并探讨如何利用模拟退火器在普通电脑上运行或真实的量子退火服务来求解。你会发现即使没有量子硬件这套建模思想和工具链本身也能为我们解决复杂优化问题打开一扇新窗。2. 问题拆解投资组合优化的数学内核与QUBO转化逻辑在动手写代码之前我们必须把问题吃透。一个标准的均值-方差投资组合优化模型其数学核心通常包含以下几个部分决策变量对于N种资产我们需要决定每个资产的权重 ( w_i )投资比例。在经典连续模型中( w_i ) 是0到1之间的实数。但为了适配QUBO其变量是二进制的0或1我们通常需要做离散化处理。一个常见技巧是使用多个二进制变量来表示一个权重。例如如果我们希望权重精度到0.01即1%那么对于每种资产i我们可以用K个二进制变量 ( x_{i, k} ) 来表示其中 ( k0,1,...,K-1 )。那么权重 ( w_i ) 可以表示为 [ w_i \frac{1}{C} \sum_{k0}^{K-1} 2^k x_{i, k} ] 这里的C是一个归一化常数确保所有权重之和为1。这种方法叫二进制编码它把连续问题转化为了离散组合问题。目标函数风险最小化马科维茨理论的核心是方差代表风险。投资组合的方差风险可以表示为 ( \mathbf{w}^T \Sigma \mathbf{w} )其中 ( \mathbf{w} ) 是权重向量( \Sigma ) 是资产的协方差矩阵。这是一个二次型天然符合QUBO的形式QUBO的目标函数就是 ( \sum_{i,j} Q_{ij} x_i x_j ) 的形式。约束条件这是将现实问题转化为无约束QUBO模型的关键也是PyQUBO大显身手的地方。主要约束通常有预算约束所有资产权重之和为1即 ( \sum_{i1}^N w_i 1 )。收益约束投资组合的期望收益不低于某个目标值R即 ( \sum_{i1}^N \mu_i w_i \geq R )其中 ( \mu_i ) 是资产i的期望收益。整数约束/边界约束比如不允许卖空( w_i \geq 0 )或者单个资产权重有上限下限。如何将约束“塞进”QUBOQUBO本身是无约束的优化形式。PyQUBO采用惩罚函数法来处理约束。基本原理是将约束条件转化为一个惩罚项加到原始目标函数上。如果解违反了约束这个惩罚项就会产生一个很大的正值从而提高目标函数值因为我们求的是最小值使得这个解变得“不优”。例如对于预算约束 ( \sum_i w_i - 1 0 )我们可以将其转化为惩罚项 ( \lambda (\sum_i w_i - 1)^2 )其中 ( \lambda ) 是一个很大的正数称为惩罚系数。当约束被满足时此项为0不满足时此项为正起到惩罚作用。PyQUBO的魅力在于你只需要用Python语法写出这个约束表达式它会自动帮你展开、化简并合并到整体的QUBO系数矩阵中。注意惩罚系数 (\lambda) 的选取是个技术活。太小了约束不起作用解会乱跑太大了可能会掩盖原始目标函数风险最小化的细节导致找到的解虽然满足约束但风险并非最优。通常需要根据目标函数的量级进行试验和调整。一个经验法则是让惩罚项的数量级显著大于目标函数可能的变化范围。所以整个建模的思维链条是现实问题 - 数学建模含约束- 二进制编码离散化 - 用PyQUBO定义目标函数和约束自动生成惩罚项- 编译为QUBO系数矩阵。接下来我们就用代码把这个链条实现出来。3. 手把手PyQUBO建模从资产数据到QUBO矩阵假设我们有5种资产的历史收益率数据我们想构建一个投资组合。首先我们需要计算它们的期望收益向量mu和协方差矩阵sigma。这里我们用随机数生成来模拟实战中你需要替换为真实数据。import numpy as np import pandas as pd from pyqubo import Array, Constraint, Placeholder # 1. 模拟资产数据 np.random.seed(42) n_assets 5 # 生成期望年化收益率 (例如在 -0.1 到 0.2 之间) mu np.random.uniform(-0.1, 0.2, n_assets) # 生成协方差矩阵 (确保是半正定的) random_matrix np.random.randn(n_assets, n_assets) sigma np.dot(random_matrix, random_matrix.T) / 10 # 缩放以控制风险量级 print(期望收益向量 mu:\n, mu) print(\n协方差矩阵 sigma:\n, sigma)接下来是核心的建模部分。我们设定用4个二进制变量K4来表示每种资产的权重精度可以达到 ( 1/(2^4 -1) \approx 0.067 )。权重 ( w_i \frac{1}{C} \sum_{k0}^{3} 2^k x_{i,k} )其中C是归一化常数为了满足权重和为1我们可以让 ( C \sum_{i0}^{4} \sum_{k0}^{3} 2^k x_{i,k} )但这会使表达式复杂。更常用的方法是先定义未归一化的权重 ( u_i \sum_{k0}^{3} 2^k x_{i,k} )然后在约束中要求 ( \sum_i u_i S )其中S是一个我们设定的表示“总份额”的常数比如 ( 2^K -1 )最后实际权重 ( w_i u_i / S )。这样建模更清晰。# 2. 定义QUBO变量和参数 K 4 # 每种资产用4个比特表示 S 2**K - 1 # 总份额当所有比特为1时u_i的最大值 # 创建二进制变量数组形状为 (n_assets, K) x Array.create(x, shape(n_assets, K), vartypeBINARY) # 定义未归一化的权重 u_i u [sum(2**k * x[i, k] for k in range(K)) for i in range(n_assets)] # 3. 构建目标函数投资组合方差 (风险) # 注意现在我们的变量是u需要转换为实际权重 w_i u_i / S # 组合方差 w^T sigma w (1/S^2) * u^T sigma u # 因为1/S^2是常数最小化 u^T sigma u 等价于最小化方差 H_risk 0 for i in range(n_assets): for j in range(n_assets): H_risk sigma[i, j] * u[i] * u[j] # 4. 构建约束条件 # 4.1 预算约束所有未归一化权重之和等于总份额S (即归一化后权重和为1) budget_constraint Constraint((sum(u) - S)**2, labelbudget) # 4.2 收益约束期望收益 目标收益R # 组合收益 sum(mu_i * w_i) (1/S) * sum(mu_i * u_i) R 0.08 # 目标年化收益8% # 注意约束是 sum(mu_i * u_i) / S R 转化为 sum(mu_i * u_i) R * S # 在QUBO中处理“大于等于”约束通常引入松弛变量或转化为等式。 # 一个技巧是定义收益惩罚项为 (R*S - sum(mu_i * u_i))^2但只有当 sum(mu_i * u_i) R*S 时才惩罚。 # 然而平方项对于“大于等于”是不对称的。更稳妥的方法是引入一个二进制松弛变量y_l将不等式转化为等式。 # 这里为了简化演示我们采用惩罚函数形式并设置一个较大的惩罚系数但这可能不是最优方法。 # 更严谨的做法参见后续的“高级约束处理”部分。此处先做简化 H_return_penalty (R * S - sum(mu[i] * u[i] for i in range(n_assets)))**2 # 我们将其作为一个惩罚项加入而不是用Constraint对象以便单独控制其系数。 # 实际上对于“”约束更好的形式是 max(0, R*S - sum(...))^2但这在QUBO中是非二次的。 # 因此对于生产环境建议使用引入辅助变量的方法。 # 5. 构建总哈密顿量 (QUBO模型) # Placeholder 用于在编译后动态设置惩罚系数 lambda_budget Placeholder(lambda_budget) lambda_return Placeholder(lambda_return) H H_risk lambda_budget * budget_constraint lambda_return * H_return_penalty # 6. 编译模型 model H.compile()到这里我们已经用PyQUBO定义好了整个问题。model对象包含了我们模型的内部表示。接下来我们需要将其“编译”成真正的QUBO系数字典也就是一个映射(i, j)到系数 ( Q_{ij} ) 的字典。同时我们需要给占位符Placeholder赋予具体的值。# 7. 生成QUBO系数 # 设定惩罚系数这是一个需要调试的参数 feed_dict {lambda_budget: 1e4, lambda_return: 1e3} qubo, offset model.to_qubo(feed_dictfeed_dict) # qubo 是一个字典键是变量索引的元组 (i, j)值是它们之间的二次项系数。 # 例如qubo[(0,1)] 是变量 x_0 和 x_1 的相互作用系数。 # offset 是一个常数项在比较能量时有用但通常不影响最优解的顺序。 print(fQUBO字典中的交互项数量: {len(qubo)}) print(f模型中的变量总数: {model.num_binary_variables}) # 注意变量总数是 n_assets * K 5 * 4 20现在qubo这个字典就是我们可以喂给各种退火求解器无论是模拟的还是量子的的“标准餐食”。offset是常数项在求解时通常可以忽略但在计算最终目标函数值时需要加上。4. 求解与解码模拟退火、量子退火与结果分析得到QUBO模型后我们就可以求解了。由于我们大多数人没有直接的量子退火机可用我们先使用经典的模拟退火算法它可以在普通计算机上运行原理是模拟固体退火过程通过缓慢降低“温度”来寻找低能态即QUBO的最小值。dimod库提供了方便的模拟退火求解器。import dimod # 1. 使用dimod的模拟退火求解器 sampler dimod.SimulatedAnnealingSampler() # 我们需要将pyqubo的qubo格式转换为dimod的BinaryQuadraticModel bqm dimod.BinaryQuadraticModel.from_qubo(qubo, offsetoffset) # 运行模拟退火可以指定读取数num_reads和退火参数 response sampler.sample(bqm, num_reads1000, beta_range[0.1, 20.0], num_sweeps1000) # 2. 查看结果 print(模拟退火结果示例能量最低的5个解:) print(response.truncate(5))response是一个包含多个样本解的对象每个样本都有对应的变量赋值和能量值即QUBO目标函数值。能量越低解越好。我们需要从中找出能量最低且满足约束的解尽管有惩罚项但有时惩罚系数不足解可能轻微违反约束我们需要检查。# 3. 解码并验证最佳解 # 获取能量最低的解 best_sample response.first.sample # 这是一个变量索引到0/1值的字典 best_energy response.first.energy # 将二进制解解码回权重 # 首先我们需要将扁平的变量字典还原成二维数组 x[i][k] x_solution np.zeros((n_assets, K), dtypeint) var_list list(model.variables) # 获取变量名列表如 [x[0][0], x[0][1], ...] # 构建一个从变量名到索引的映射这里索引是我们在qubo中用的数字索引 # 注意model.variables 的顺序可能与 to_qubo 生成的索引顺序一致。 # 更稳妥的方法是利用 model.decode_sample best_solution_decoded, broken_constraints, energy model.decode_sample(best_sample, vartypeBINARY, feed_dictfeed_dict) print(f\n解码后的信息:) print(f 能量 (包含惩罚项): {energy}) print(f 违反的约束: {broken_constraints}) # 我们希望这个是空的 # 提取未归一化的权重 u u_solution [sum(2**k * best_solution_decoded[x][i][k] for k in range(K)) for i in range(n_assets)] # 计算实际权重 w w_solution np.array(u_solution) / S # 计算实际组合收益和风险 portfolio_return np.dot(mu, w_solution) portfolio_risk np.dot(w_solution, np.dot(sigma, w_solution)) print(f\n最优投资组合权重:) for i in range(n_assets): print(f 资产 {i}: {w_solution[i]:.4f} ({w_solution[i]*100:.2f}%)) print(f 权重总和: {np.sum(w_solution):.6f} (应接近1)) print(f\n组合预期年化收益: {portfolio_return:.4f} ({portfolio_return*100:.2f}%)) print(f组合预期年化风险 (方差): {portfolio_risk:.6f}) print(f目标收益R: {R:.4f})运行这段代码你就能看到一个由模拟退火找到的投资组合分配方案。如果broken_constraints不为空或者权重和严重偏离1或收益低于目标你可能需要回头调整惩罚系数lambda_budget和lambda_return然后重新编译QUBO并求解。那么量子退火呢如果你有D-Wave的访问权限过程非常相似。你需要安装dwave-ocean-sdk配置好API令牌和求解器端点然后将bqm提交给量子退火器。# 示例使用D-Wave量子退火需要配置环境 from dwave.system import DWaveSampler, EmbeddingComposite # sampler_q EmbeddingComposite(DWaveSampler(solver{qpu: True})) # response_q sampler_q.sample(bqm, num_reads1000) # 后续解码分析过程与模拟退火完全相同量子退火器的优势在于对于某些特定结构的QUBO问题它可能通过量子隧穿效应更有效地跳出局部最优解从而有可能找到比经典模拟退火更好的解。但对于小规模问题经典算法通常已经足够。5. 高级话题与实战避坑指南走到这一步你已经能用PyQUBO解决一个基本的投资组合优化问题了。但在实际应用或比赛中你会遇到更多挑战。下面分享几个关键的高级话题和踩坑经验。5.1 不等式约束的精确处理前面我们用简单的平方惩罚来处理收益约束sum(mu_i * u_i) R*S这并不精确因为它对超过目标收益的情况也进行了惩罚虽然惩罚较小。标准的处理方法是引入辅助变量Slack Variable。例如将不等式A B改写为等式A - B - s 0其中s是一个非负的整数或连续变量。为了放入QUBO我们需要对s进行二进制编码。假设我们预估A-B的最大值为M我们可以用L个二进制位来表示ss sum_{l0}^{L-1} 2^l * y_l。这样约束A - B - s 0就可以用(A - B - s)^2作为惩罚项精确地表达了。在PyQUBO中你需要额外创建y_l这些二进制变量。这增加了变量总数但保证了约束的精确性。5.2 二进制编码的精度与规模权衡我们用K个比特表示一种资产的权重这决定了权重的分辨率最小变化单位是 (1/S)。K越大精度越高但QUBO问题的变量数呈线性增长N * K问题规模迅速膨胀。量子退火机如D-Wave的量子比特数是有限的且并非所有比特之间都能直接耦合。变量越多问题越难嵌入到物理硬件中求解也越困难。实战建议从较小的K如3或4开始验证模型逻辑。根据资产数量和可用计算资源逐步增加K。对于大规模资产N50可能需要考虑更高级的编码方式如对数编码或先进行资产聚类降维。5.3 惩罚系数调参的艺术惩罚系数lambda的选择极大影响求解效果。我的经验是分层调试先只加一个约束如预算约束调整其lambda直到该约束被完全满足broken_constraints为空。记录下这个lambda的大致量级。叠加调试加入第二个约束如收益约束将其lambda设为一个较小的值比如第一个约束lambda的1/10然后逐步增大观察两个约束是否同时被满足以及目标函数值风险的变化。观察能量尺度检查最终QUBO中各项系数的数量级。理想情况下惩罚项产生的系数应显著大于目标函数风险项的系数但又不能大到让风险项的差异被完全淹没。通常惩罚项系数比风险项大1到3个数量级是一个不错的起点。利用PlaceholderPyQUBO的Placeholder功能就是为了方便调参设计的。你可以编译一次模型然后快速生成不同惩罚系数下的QUBO进行测试无需重新编译整个模型。5.4 从QUBO解到实际权重的后处理由于二进制编码和离散化我们得到的最优解对应的权重w_i可能不是数学上绝对的最优解因为连续空间被离散化了。一个常见的后处理技巧是固定大部分资产的权重只对少数权重在离散点附近进行微调。例如你可以将QUBO解作为初始解在其附近用一个快速的连续优化器如SciPy的minimize对原始连续问题进行局部搜索但需要满足离散解已满足的约束。这通常能进一步提升解的质量。5.5 量子退火的实际考量如果你真的使用D-Wave嵌入问题不是所有QUBO都能直接映射到硬件的耦合图上。D-Wave的Ocean SDK提供的EmbeddingComposite会自动处理嵌入但这会引入额外的辅助量子比特并可能影响性能。对于大规模问题可能需要手动设计嵌入或使用分治算法。读取数与退火参数num_reads相当于从退火机读取多少个样本。越多找到好解的概率越大但也越耗时费钱D-Wave按QPU时间计费。annealing_time、chain_strength等参数需要根据具体问题调整Ocean SDK有工具帮助自动调参。退火与反向退火标准退火是从初始状态开始。反向退火则允许你从一个经典解比如用模拟退火找到的好解开始让量子退火机在其周围进行局部搜索这常用于解决方案的精细化。通过这个完整的流程——从问题理解、数学建模、PyQUBO编码、模拟/量子求解到结果分析与高级调优——你不仅掌握了解决一个具体赛题的工具更获得了一套将经典组合优化问题转化为QUBO模型并利用前沿计算平台求解的通用方法论。这其中的思维转换和工程实践细节才是应对未来更多复杂优化挑战的真正武器。