1. 项目概述当量子计算遇上金融风控去年带队参加MathorCup选的就是这道A题。说实话当时看到“量子计算机”和“信用评分卡”这两个词放在一起团队里的小伙伴第一反应都是懵的。一个听起来是前沿物理实验室的玩意儿另一个是银行信贷部门天天打交道的传统工具这俩能扯上什么关系但正是这种跨领域的碰撞让这道题充满了挑战和魅力。它本质上不是让你去造一台量子计算机而是要求你理解一种名为“量子退火”的专用计算范式并学会如何将一个经典的、复杂的金融优化问题翻译成这种新硬件能“听懂”的语言——也就是QUBO模型。这个过程就像是为一个只会意大利语的歌剧演唱家量子退火器准备一份他能完美演绎的乐谱QUBO模型而乐谱的内容则是如何从成千上万张信用卡申请中科学地挑选出那些既能控制风险又能最大化利润的客户。这道题的核心目标非常明确在给定的风险资本约束下从多张信用评分卡中选出一个最优的卡片组合并为每张选中的卡设定一个最优的通过率阈值最终实现整体利润的最大化。这实际上是一个经典的组合优化问题但变量多、约束复杂用传统算法求解可能面临“组合爆炸”的困境。量子计算特别是量子退火因其在解决某些特定组合优化问题上的潜在优势而被引入作为探索方案。整个项目的价值不仅在于求解一道赛题更在于提供了一个绝佳的实践窗口让我们这些建模爱好者能亲手触摸到“量子计算金融科技”这个前沿交叉领域的门槛理解从业务问题抽象到数学模型再适配到新型计算架构的全链条思维。2. 问题拆解从金融业务到数学语言面对一个复杂问题直接上手编码是最大的忌讳。我们必须先抽丝剥茧把模糊的业务描述转化成精确的数学定义。这是建模过程中最考验基本功也最决定成败的一步。2.1 核心要素定义首先我们需要明确题目给出的所有“积木”信用评分卡假设有N张。每张卡i可以独立地对客户进行信用评分本质上是一个分类器。通过率阈值对于每张卡i我们可以设定一个阈值theta_i。评分高于此阈值的客户通过审批低于则拒绝。theta_i是一个连续变量通常在0到1之间表示分数分位点。坏账率这是核心风险指标。对于卡i在给定阈值theta_i下通过审批的客户中最终违约的比例记为b_i(theta_i)。通常阈值设得越高通过率越低通过的客户质量越好坏账率b_i越低。这是一个关于theta_i的单调递减函数。通过人数在阈值theta_i下通过审批的客户数量记为g_i(theta_i)。这是一个关于theta_i的单调递减函数。利润每通过一个客户如果其正常还款银行能获得固定收益r如果其违约银行会损失本金L。因此对于卡i在阈值theta_i下单张卡的总利润P_i(theta_i)可以表示为P_i(theta_i) g_i(theta_i) * [ r*(1 - b_i(theta_i)) - L*b_i(theta_i) ]风险资本银行需要为可能发生的坏账计提资本。题目约束总风险资本不超过C。通常风险资本与坏账总额成正比即Risk_i(theta_i) g_i(theta_i) * L * b_i(theta_i) * k其中k是资本计提系数。总风险资本为各卡风险资本之和。2.2 决策变量与问题形式化我们的决策有两层选择哪些卡用一个二进制变量x_i表示x_i 1表示选择第i张卡x_i 0表示不选。每张选中的卡设定什么阈值连续变量theta_i。那么总利润和总风险资本分别为总利润 P_total sum_{i1}^{N} [ x_i * P_i(theta_i) ]总风险资本 R_total sum_{i1}^{N} [ x_i * Risk_i(theta_i) ]优化问题可以表述为最大化P_total 约束条件R_total C 决策变量x_i ∈ {0, 1}, theta_i ∈ [0, 1] (且当 x_i0 时theta_i 无意义可强制为0或忽略)这是一个混合整数非线性规划问题。x_i是整数theta_i是连续变量目标函数和约束条件由于g_i(theta_i)和b_i(theta_i)的存在而高度非线性。注意这里有一个关键的建模技巧。theta_i是连续的但量子退火器处理离散变量更拿手。因此在实际操作中我们几乎总是需要对theta_i进行离散化。例如将[0,1]区间均匀划分为M个等级如0.05, 0.10, ..., 1.00这样theta_i就变成了一个离散选择可以用一组辅助的二进制变量来表示。这是将问题转化为QUBO的关键一步。2.3 引入QUBO模型框架量子退火器如D-Wave原生求解的是二次无约束二进制优化问题其标准形式如下最小化y sum_{i} a_i q_i sum_{ij} b_{ij} q_i q_j 其中q_i ∈ {0, 1}我们的任务就是把那个带有约束的混合整数非线性规划问题变形到这个简单的二次多项式形式上。这需要两大步骤离散化和约束处理。离散化连续变量对于每张卡i和每个离散化的阈值选项m(共M个)我们引入一个二进制变量q_{i,m}。q_{i,m} 1表示“为第 i 张卡选择第 m 个阈值”。为了保证每张卡至多选择一个阈值或不选即所有q_{i,m}0我们需要添加惩罚项A * (sum_{m} q_{i,m} - x_i)^2其中A是一个很大的正数惩罚系数。当且仅当选择的阈值数等于x_i0或1时这个惩罚项为0。处理不等式约束风险资本约束R_total C是一个不等式。我们将其转化为等式R_total s C其中s 0是松弛变量。松弛变量也需要被离散化并用二进制变量表示。然后将等式约束以惩罚项形式加入目标函数B * (R_total s - C)^2其中B是另一个很大的惩罚系数。这样当约束被违反时目标函数值会急剧增大迫使优化器寻找满足约束的解。最终我们的最大化总利润问题转化为了一个最小化包含原始利润项取负和一系列惩罚项的二次型问题。这就是完整的QUBO模型。3. 数据预处理与函数拟合在实际比赛或应用中我们通常不会直接得到g_i(theta)和b_i(theta)的解析式而是拥有每个评分卡在历史样本上的数据每个客户的评分和最终是否违约的标签。3.1 构建阈值-性能曲线我们需要从数据中提炼出这两个关键函数。步骤如下排序对于每张卡i将其所有历史客户按评分从高到低排序。滑动阈值从最高分最严格的通过率例如前1%开始逐步滑动到最低分通过率100%。计算指标在每个阈值点theta对应一个通过率计算通过人数 g_i(theta)评分 当前阈值的客户数。坏账率 b_i(theta)在通过客户中违约客户的数量 / 通过总人数。数据存储得到一系列离散的(theta, g_i, b_i)数据点。3.2 函数拟合与插值为了在优化中能够计算任意theta下的g_i和b_i我们需要对离散数据点进行拟合或插值。拟合可以尝试用指数函数、对数函数或分段线性函数去拟合b_i(theta)曲线通常单调递减。g_i(theta)本质上是通过率的线性函数如果总客户数固定相对简单。插值更稳健的方法是使用样条插值。例如使用三次样条插值它能保证曲线的平滑性这对于后续优化求解的稳定性很重要。Python的scipy.interpolate库中的CubicSpline或UnivariateSpline非常适用于此。import numpy as np from scipy.interpolate import CubicSpline # 假设 thresholds 是离散化的通过率阈值点如 [0.01, 0.02, ..., 1.0] # bad_rates 是对应的坏账率数组 # total_clients 是该评分卡历史客户总数 thresholds np.array([...]) bad_rates np.array([...]) passing_clients total_clients * thresholds # 近似线性关系 # 创建坏账率关于阈值的样条插值函数 bad_rate_spline CubicSpline(thresholds, bad_rates, bc_typenatural) # ‘natural’ 样条边界条件 # 现在可以计算任意阈值 t 下的坏账率 t 0.35 estimated_bad_rate bad_rate_spline(t) estimated_passing total_clients * t实操心得拟合时一定要检查曲线的单调性。理论上坏账率应随阈值升高而降低。如果拟合后的曲线在某些区间出现“反弹”说明拟合函数选择不当或数据有噪声这时分段线性插值可能是更安全的选择。此外对于通过人数g_i(theta)如果总客户数很大直接使用total_clients * theta是合理的简化无需复杂拟合。4. QUBO模型构建详解这是整个项目的核心技术环节。我们将一步步把业务问题“编码”成QUBO矩阵。4.1 变量映射假设有N3张卡每张卡的阈值离散化为M5个选项例如通过率20%40%60%80%100%。那么我们需要N * M 15个变量来表示阈值选择q_{0,0}, q_{0,1}, ..., q_{2,4}。还需要N3个变量来表示卡片是否被选择吗实际上不需要。卡片被选择等价于为该卡选择了某一个阈值。我们可以定义如果一张卡的所有q_{i,m}都为0则表示该卡未被选中。但这样无法直接表达“不选卡”带来的收益为0。在目标函数中不选卡自然就是利润贡献为0。为了处理风险资本约束R_total C我们引入一个离散化的松弛变量s。假设风险资本最大值为C_max我们将[0, C_max]区间离散化为K份用K个二进制变量s_k的线性组合来表示s例如采用二进制编码或one-hot编码。为简化这里假设使用one-hot编码即s sum_{k} (value_k * s_k)且sum_{k} s_k 1。因此总二进制变量数量为N*M K。4.2 目标函数构建利润项对于每个q_{i,m}如果它被激活1则带来利润P_i(theta_m)。由于QUBO是最小化问题我们将最大化利润转化为最小化负利润。因此目标函数中的线性项系数为a_{i,m} -P_i(theta_m)这里没有直接的二次项。所以利润部分在QUBO矩阵中体现为对角线上的元素。4.3 约束条件转化为惩罚项这是构建QUBO最精妙也最需要小心的地方。惩罚系数A和B的选择至关重要。1. 单卡单阈值约束对于每张卡i至多只能选择一个阈值。这可以表述为sum_{m} q_{i,m} 1。在QUBO中我们将其转化为惩罚项Penalty_single A * sum_{i} (sum_{m} q_{i,m} - 1)^2注意这里我们允许sum_{m} q_{i,m} 0不选该卡但强制其不能大于1。展开这个平方项会产生q_{i,m}的线性项对角线和q_{i,m} * q_{i,n}(m!n) 的二次项非对角线。2. 风险资本约束R_total C。首先计算选择特定阈值带来的风险资本R_i(theta_m)。定义R_{i,m} Risk_i(theta_m)。 总风险资本R_total sum_{i} sum_{m} (R_{i,m} * q_{i,m})。 引入松弛变量s用one-hot编码的变量s_k表示其取值为value_k约束变为sum_{i,m} (R_{i,m} * q_{i,m}) sum_{k} (value_k * s_k) C将其作为惩罚项加入Penalty_risk B * ( sum_{i,m} (R_{i,m} * q_{i,m}) sum_{k} (value_k * s_k) - C )^2展开这个平方项将产生涉及q和s的所有变量的线性项和二次项。3. 松弛变量单激活约束如果使用one-hot编码需要保证sum_{k} s_k 1。这同样是一个惩罚项Penalty_slack C * (sum_{k} s_k - 1)^24.4 组装QUBO矩阵最终我们的QUBO问题形式为最小化H (-利润项) Penalty_single Penalty_risk Penalty_slack将上述所有项展开、合并同类项我们就能得到一个对称的或上三角QUBO矩阵Q其中Q[j][j]是变量j的线性系数Q[j][k](jk) 是变量j和k的二次项系数。优化目标就是找到二进制向量q使得q^T Q q最小。import numpy as np def build_qubo_matrix(N, M, K, profits, risks, slack_values, C, A, B, C_slack): 构建QUBO矩阵 N: 卡片数量 M: 每张卡的阈值选项数 K: 松弛变量选项数 profits: 形状为 (N, M) 的数组profits[i][m] 表示卡i选择阈值m的利润 risks: 形状为 (N, M) 的数组risks[i][m] 表示对应的风险资本 slack_values: 长度为K的数组松弛变量各选项的值 C: 风险资本上限 A, B, C_slack: 惩罚系数 total_vars N * M K Q np.zeros((total_vars, total_vars)) # 1. 负利润项线性项放在对角线 for i in range(N): for m in range(M): idx i * M m Q[idx][idx] -profits[i][m] # 最小化负利润 最大化利润 # 2. 单卡单阈值约束惩罚项展开 # 对于每张卡i惩罚项 A * (sum_m q_{i,m} - 1)^2 A*(sum_m q_{i,m}^2 2*sum_{mn} q_{i,m}q_{i,n} - 2*sum_m q_{i,m} 1) # 常数1可以忽略。q_{i,m}^2 q_{i,m} (因为二进制变量) for i in range(N): start_idx i * M # 线性部分A * (1 - 2*1) * q_{i,m}? 让我们仔细展开 # (sum_m q_{i,m} - 1)^2 sum_m q_{i,m}^2 2*sum_{mn} q_{i,m}q_{i,n} - 2*sum_m q_{i,m} 1 # 对于每个固定的i和m # q_{i,m}^2 的系数是 A加到对角线 Q[idx][idx] # q_{i,m} 的系数是 -2A加到对角线 Q[idx][idx] (因为线性项在对角线) # q_{i,m} * q_{i,n} 的系数是 2A加到 Q[idx_m][idx_n] (m n) for m in range(M): idx_m start_idx m Q[idx_m][idx_m] A * (1 - 2) # 来自 q^2 和 -2q 项 for n in range(m1, M): idx_n start_idx n Q[idx_m][idx_n] 2 * A # 3. 风险资本约束惩罚项展开 # 惩罚项 B * (sum_{i,m} R_{i,m}*q_{i,m} sum_k val_k*s_k - C)^2 # 令 S sum_{i,m} R_{i,m}*q_{i,m} sum_k val_k*s_k - C # 则 B * S^2 B * [ (第一部分)^2 (第二部分)^2 2*(第一部分)*(第二部分) C^2项 - 2C*(第一部分第二部分) ] # 常数项C^2可忽略。我们分别处理。 # 先处理涉及q变量的部分 for i1 in range(N): for m1 in range(M): idx1 i1 * M m1 R1 risks[i1][m1] # 线性项来自 -2C*R1 和 B*R1^2 (来自平方项但q是二进制q^2q) Q[idx1][idx1] B * (R1*R1 - 2*C*R1) # 二次项不同q变量之间的相互作用 2B * R1 * R2 for i2 in range(N): start_m2 0 if i2 i1 else m11 # 避免重复计算 for m2 in range(start_m2, M): if i2 i1 and m2 m1: continue idx2 i2 * M m2 R2 risks[i2][m2] Q[idx1][idx2] 2 * B * R1 * R2 # 处理松弛变量s部分 (索引从 N*M 开始) s_start N * M for k1 in range(K): idx1 s_start k1 val1 slack_values[k1] # 线性项 Q[idx1][idx1] B * (val1*val1 - 2*C*val1) # 二次项s变量之间 for k2 in range(k11, K): idx2 s_start k2 val2 slack_values[k2] Q[idx1][idx2] 2 * B * val1 * val2 # 二次项q变量与s变量之间 for i in range(N): for m in range(M): idx_q i * M m R risks[i][m] Q[idx_q][idx1] 2 * B * R * val1 # 注意索引顺序保证 idx_q idx1 # 4. 松弛变量单激活约束 for k1 in range(K): idx1 s_start k1 # 线性项来自 C_slack * (1 - 2*1) * s_k Q[idx1][idx1] C_slack * (1 - 2) for k2 in range(k11, K): idx2 s_start k2 Q[idx1][idx2] 2 * C_slack return Q关键技巧惩罚系数A, B, C_slack的选取这是QUBO建模的“艺术”部分。系数太小约束得不到严格执行系数太大可能掩盖原始目标利润导致求解器只满足约束而不管利润高低。经验法则惩罚系数应显著大于目标函数中变量的典型取值范围。例如利润项的数量级在1e4左右那么惩罚系数可以从1e5或1e6开始尝试。分层设置通常A单选择约束可以设得最大因为它是最硬的约束。B风险约束次之。C_slack松弛变量约束可以设得和A类似或稍小。测试验证必须用一些小规模实例进行测试。随机生成一些解手动计算总惩罚项和原始目标确保违反约束的解对应的QUBO能量值远高于可行解的能量值。一个常用的方法是系数 |最大可能利润差距| / |最小约束违反度|。5. 求解与后处理从量子退火到经典求解5.1 量子退火求解流程如果你有权限访问真实的量子退火器如通过D-Wave的Leap云服务流程如下问题嵌入QUBO矩阵需要映射到量子退火器的物理量子比特连接图Chimera或Pegasus拓扑上。这个过程称为“嵌入”可以使用D-Wave提供的minorminer或EmbeddingComposite工具自动完成。对于完全连接图我们的QUBO矩阵通常是稠密的需要多个物理量子比特链式连接来代表一个逻辑变量这会消耗大量资源。参数调优设置退火时间、读取次数等参数。较长的退火时间通常有助于找到更优解但会增加计算成本。采样求解提交任务并获取多个解样本。解译将返回的二进制样本向量映射回我们定义的变量q_{i,m}和s_k。# 伪代码示例使用D-Wave Ocean SDK from dwave.system import DWaveSampler, EmbeddingComposite import dimod # 假设我们已经得到了QUBO矩阵 Q (作为上三角字典或dimod.BQM对象) bqm dimod.BinaryQuadraticModel.from_qubo(Q) # 将矩阵转换为BQM对象 # 使用嵌入复合器处理硬件连接问题 sampler EmbeddingComposite(DWaveSampler(tokenYOUR_TOKEN, solverAdvantage_system6.4)) # 执行退火读取多个样本 sampleset sampler.sample(bqm, num_reads1000, annealing_time100) # 查看最优解 best_sample sampleset.first.sample best_energy sampleset.first.energy5.2 经典模拟退火求解对于大多数参赛者使用真实的量子硬件并不现实。幸运的是我们可以使用经典模拟退火算法来求解QUBO模型其原理类似且有很多成熟的库。import neal # 使用neal库的模拟退火求解器 sampler neal.SimulatedAnnealingSampler() sampleset sampler.sample(bqm, num_reads1000, num_sweeps1000) best_sample sampleset.first.sample模拟退火的结果可以作为基准用于验证后续量子求解的结果或者直接作为最终解决方案。5.3 解的后处理与验证获得最优的二进制向量best_sample后我们需要将其解码为业务决策解析卡片选择与阈值对于每张卡i检查q_{i,0}到q_{i, M-1}。有且仅有一个为1则选中该卡并记录其对应的阈值theta_m。如果全部为0则该卡未被选中。解析松弛变量检查s_k为1的那个k对应的slack_values[k]就是松弛变量s的值。验证约束重新计算选中卡的总风险资本R_total加上松弛变量s验证是否等于在数值误差内约束上限C。计算最终利润根据选中的卡和阈值使用原始数据或拟合函数计算总利润P_total。注意事项由于退火过程的随机性我们可能会得到多个能量相近的解。一个好的实践是分析前几十个最优解观察它们结构上的共性例如是否某几张卡总是被选中这能增加方案的可信度。同时一定要用经典优化器如Gurobi, CPLEX针对简化后的MIP模型或暴力枚举对于小规模问题进行交叉验证确保QUBO模型和求解过程没有根本性错误。6. 模型拓展与优化思考在实际比赛中要脱颖而出还需要在基础模型上做深度思考和拓展。6.1 非线性利润与风险函数的处理我们之前假设g_i(theta)是线性的b_i(theta)通过插值得到。但如果数据表现出强烈的非线性或者利润函数本身不是简单的[r*(1-b) - L*b]而是更复杂的公式例如包含资金成本、运营成本我们的模型依然可以容纳。只需要在计算P_i(theta_m)和R_i(theta_m)时使用更复杂的函数即可。QUBO模型本身不关心这些值是如何来的它只关心每个决策变量组合对应的最终数值。6.2 多周期动态优化题目可以拓展为多周期例如12个月的优化。此时决策变量可能包含每个周期是否使用某张卡、以及阈值如何随时间调整。约束可能包括风险资本的跨周期分配、客户群体的动态变化等。这会将问题复杂度提升一个数量级但建模框架不变定义时间索引t变量变为q_{i,m,t}目标函数变为多期利润的净现值总和约束也需要按时间展开。QUBO模型依然适用但变量数会剧增对求解器规模要求更高。6.3 与传统优化方法的对比分析在论文中一个重要的部分是对比。除了量子模拟退火这个问题显然可以用传统的数学规划方法求解例如混合整数非线性规划使用IPOPT、Bonmin等求解器。转化后的混合整数线性规划如果对g_i和b_i进行分段线性化可以将问题近似为MILP问题用Gurobi、CPLEX高效求解。元启发式算法如遗传算法、粒子群算法直接优化x_i和theta_i。你需要设计实验在相同数据和小规模场景下对比这些方法在求解质量最优利润、求解速度和可扩展性变量增多时的表现上的差异。量子退火可能在中等规模、特定结构的问题上展现出潜力但在小问题上可能不如传统MILP精确在大问题上受限于硬件量子比特数和连接性。客观地分析优劣是论文的亮点。6.4 实际部署的考量如果真的在银行系统中应用还需要考虑冷启动问题新上线的评分卡没有历史数据如何估计b_i(theta)可能需要基于业务经验设定保守参数或利用迁移学习从类似卡片获取先验知识。在线学习与动态调整模型应该定期如每月用最新数据重新训练和优化以适应客户群体和宏观经济环境的变化。解释性量子退火模型是一个“黑箱”如何向风控委员会解释为什么选择这组卡片和阈值需要结合SHAP等可解释性AI工具或者准备一个基于规则的、近似最优的备用方案。7. 代码实现要点与避坑指南7.1 完整代码结构一个健壮的实现应该包含以下模块data_loader.py: 读取历史评分卡数据。curve_fitting.py: 拟合或插值得到每张卡的g_i(theta)和b_i(theta)函数。precompute.py: 根据离散化的阈值网格预计算所有P_i(theta_m)和R_i(theta_m)。qubo_builder.py: 构建QUBO矩阵的核心函数包含惩罚项系数的设置。solver.py: 封装模拟退火或量子退火求解器的调用。post_processor.py: 解析求解结果计算最终指标并输出可读的报告。main.py: 主流程控制器。config.yaml: 配置文件集中管理卡片数量、离散化粒度、惩罚系数、求解器参数等。7.2 常见陷阱与调试技巧陷阱一惩罚系数失衡。症状求解器返回的解总是违反约束或者总是选择很少的卡利润目标被压制。调试打印出几个手动构造的可行解和不可行解对应的QUBO能量值q^T Q q。确保可行解的能量值远低于不可行解。逐步增大惩罚系数直到满足这个条件。陷阱二离散化粒度不足。症状求解结果不稳定细微调整参数后最优解变化很大。调试增加阈值离散化的点数M。观察最优利润随M增加的变化当利润趋于稳定时说明粒度足够。陷阱三数值溢出/精度问题。症状QUBO矩阵元素过大或过小导致求解器数值不稳定。调试对利润和风险数据进行标准化。例如将利润除以一个基准值如最大单卡利润将风险资本除以约束C。这能确保QUBO矩阵中各项数量级相近。记得在最终解释结果时反标准化。陷阱四模拟退火陷入局部最优。症状多次运行得到的结果差异很大。调试增加num_reads读取次数和num_sweeps扫描次数。尝试不同的初始温度表和降温策略。使用更高级的算法如并行回火。陷阱五结果不可行。即使惩罚系数很大由于求解器的随机性返回的最优样本也可能轻微违反约束。后处理实现一个简单的修复启发式。例如如果总风险资本略微超过C可以尝试微调阈值选择相邻的更严格的阈值来降低风险直到满足约束并计算利润损失。这是一个实用的工程化步骤。7.3 性能优化建议向量化计算在构建QUBO矩阵时避免使用多层嵌套循环尽量使用NumPy的广播和向量化操作可以带来百倍的性能提升。稀疏矩阵表示对于大规模问题QUBO矩阵可能非常稀疏虽然完全连接图导致二次项稠密。使用scipy.sparse格式存储可以节省大量内存。并行计算预计算P_i和R_i对不同卡和阈值是独立的可以并行。模拟退火的多次读取 (num_reads) 也是天然并行的。完成整个项目后最大的收获不是学会了某个特定工具而是掌握了一套将复杂的、有约束的现实世界问题转化为一种特定计算模型QUBO的通用思维框架。这种“翻译”能力正是连接量子计算潜力与真实应用场景的关键桥梁。在测试中我发现对于变量数在50-100左右的问题经典模拟退火已经能在可接受时间内给出优质解而当变量规模继续增大量子退火硬件在原理上的并行性优势才会更值得期待。因此在现阶段将其视为一个强大的专用优化协处理器与传统优化方法协同工作是更务实的定位。