FPGA实现离散模拟分岔算法优化组合问题求解
1. 项目概述：FPGA实现的离散模拟分岔算法架构
在资源分配、物流调度等实际场景中，组合优化问题（Combinatorial Optimization, CO）的求解往往面临NP难问题的指数级复杂度挑战。传统CPU在处理这类问题时，随着问题规模扩大，计算时间会呈爆炸式增长。这就像要在所有可能的路线中找到最短路径，当城市数量超过30个时，传统计算机可能需要数年才能完成计算。
模拟分岔机器（Simulated Bifurcation Machines, SBMs）作为一种新型硬件加速器，通过模拟非线性参量振荡器网络的分岔现象，为组合优化问题提供了高效的解决方案。其核心优势在于算法的高度并行性，特别适合FPGA等可编程硬件实现。离散模拟分岔（discrete Simulated Bifurcation, dSB）算法作为最新变体，通过离散化处理进一步减少了模拟误差，同时保持了良好的硬件友好性。
我们开发的这个开源硬件架构，采用SystemVerilog语言描述，主要特点包括：
- 支持dSB算法及其加热机制变体
- 三级并行化设计（矩阵运算的行/列并行及模块复制）
- 可配置的并行度参数（Pc/Pr/Pb）
- 固定点数表示法优化资源利用
- 兼容max-cut和背包问题等通用伊辛模型
在AMD Kria KV260这款入门级FPGA上，我们成功实现了256变量规模的实时求解。这个开发板虽然属于低端型号（搭载Zynq UltraScale+ MPSoC），但通过架构优化，其性能已能满足许多实际应用需求。
提示：选择dSB算法而非aSB或bSB变体，主要基于三点考量：1) 离散化处理消除了模拟误差；2) 省去了乘法器资源（直接使用符号位）；3) 与加热机制兼容。这使得在同等硬件资源下，可以实现更高程度的并行化。
2. 技术背景与算法原理
2.1 伊辛模型与组合优化
伊辛模型最初用于描述磁性材料中的自旋相互作用，其哈密顿量表示为：
$$ H(s) = -\frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N J_{ij}s_is_j - B\sum_{i=1}^N h_i s_i $$
其中：
- $s_i \in {-1, +1}$ 表示第i个自旋状态
- $J_{ij}$ 为自旋间相互作用矩阵
- $h_i$ 为外部磁场影响
这个模型与组合优化中的二次无约束二进制优化（QUBO）问题等价，两者可通过简单线性变换相互转换。例如在max-cut问题中，将图节点划分为两个子集的操作，正好对应自旋的±1状态。
2.2 模拟分岔算法演进
模拟分岔算法的发展经历了三个阶段：
绝热模拟分岔（aSB）
模拟非线性Kerr振荡器的绝热演化过程，通过缓慢增加泵浦信号$a(t)$使系统发生分岔。其哈密顿量包含四次项： $$ H_{SB} = \sum_{i=1}^{n_{spin}} \left[ \frac{\Delta}{2}y_i^2 + \frac{K}{4}x_i^4 + \frac{\Delta-a(t)}{2}x_i^2 \right] - \frac{c_0}{2}\sum_{i=1}^{n_{spin}}\sum_{j=1}^{n_{spin}} J_{ij}x_i x_j $$弹道模拟分岔（bSB）
引入$x_i=\pm1$处的完全非弹性墙，消除四次项的同时减少模拟误差。当振荡器位置超出±1范围时，强制将其设置为最近的边界值。离散模拟分岔（dSB）
在bSB基础上进一步离散化，用符号函数替代位置变量： $$ f_i = \sum_{j=1}^{n_{spin}} J_{ij} \cdot \text{sgn}(x_j) $$ 这种处理虽然违反能量守恒，但能有效逃离局部极小值。
图1展示了三种算法在2节点max-cut问题中的轨迹差异。可以看到dSB的离散特性使其能快速收敛到最优解附近。
2.3 加热机制与参数选择
加热机制通过在动量更新中加入随机扰动项$\gamma y_i(t_n)\Delta t$，帮助系统跳出局部最优。这类似于模拟退火中的温度参数，但实现成本更低。
关键参数的选择策略：
- $c_0$：根据J矩阵的最大特征值$\lambda_{MAX}$设定，通常取$c_0 = \Delta/\lambda_{MAX}$
- $\Delta t$：dSB建议取1.0，bSB取0.5
- $n_{steps}$：在精度和速度间权衡，一般$10^4$步可取得较好效果
- $\gamma$：加热系数，过大导致收敛慢，过小则无法逃离局部最优
表1比较了不同算法变体在100物品背包问题上的表现：
| 算法类型 | 固定点位数 | 达到最优解概率 |
|---|---|---|
| dSB | 10-bit | 78% |
| bSB | 10-bit | 65% |
| HbSB | 15-bit | 82% |
3. 硬件架构设计
3.1 整体架构
我们的设计采用模块化结构（图2），核心组件包括：
并行计算单元（MMTE）
每个MMTE模块包含：- 矩阵-向量乘法器（MM）：处理$J_{ij}\cdot\text{sgn}(x_j)$计算
- 时间演化单元（TE）：更新$x_i$和$y_i$状态
- 本地存储器：存储当前块的自旋状态
全局存储器
- XMEM/YMEM：存储所有振荡器的位置和动量
- SGNXMEM：存储$\text{sgn}(x_i)$的符号位
- J-MEM：分块存储相互作用矩阵
控制逻辑
- 线性更新器：控制泵浦信号$a(t)$的渐变
- 调度器：协调MM和TE的流水线执行
3.2 三级并行化策略
为突破冯·诺依曼架构的串行瓶颈，我们设计了三个层次的并行：
列级并行（Pc）
单行计算中，同时处理Pc个矩阵元素。例如Pc=4时，每个周期可完成： $$ \text{acc}i += J{i,j}\text{sgn}(x_j) + J_{i,j+1}\text{sgn}(x_{j+1}) + J_{i,j+2}\text{sgn}(x_{j+2}) + J_{i,j+3}\text{sgn}(x_{j+3}) $$行级并行（Pr）
同时计算Pr行的矩阵乘法。需要复制Pr个乘法累加器（MAC），每个负责一行计算。模块级并行（Pb）
将整个系统划分为Pb个独立子块，每个子块由单独的MMTE模块处理。各模块通过共享总线同步状态。
这种设计的优势在于：
- 计算复杂度从$O(n_{spin}^2)$降至$O(n_{spin}^2/(Pc \cdot Pr \cdot Pb))$
- 资源消耗线性增长，而性能呈多项式提升
- 可根据目标FPGA资源灵活调整并行度
3.3 流水线优化
通过分析算法流程，我们发现MM和TE阶段存在天然的流水线机会：
时钟周期 | 1 | 2 | 3 | 4 | 5 | ... MM | 行1计算 | 行2计算 | 行3计算 | ... TE | 行1更新 | 行2更新 | ...要实现完美流水线，需满足： $$ \frac{n_{spin}}{Pc} \geq Pr $$ 这保证了MM计算下一批行的时间，足够TE完成当前行的更新。在256-spin系统中，选择Pr=16和Pc=16可满足该条件。
3.4 固定点数优化
软件模拟表明，10位固定点数表示已能满足dSB的精度需求（图3）。这带来两大优势：
资源节省
相比32位浮点，10位定点数：- DSP块使用减少75%
- 存储器带宽需求降低68%
- 加法器延迟缩短40%
频率提升
简化后的数据路径使时序更宽松。在KV260上，固定点设计可稳定运行在250MHz，而浮点版本仅达150MHz。
具体位宽分配：
- 整数部分：3位（覆盖[-4,4]范围）
- 小数部分：7位（精度1/128≈0.008）
4. 实现细节与优化
4.1 矩阵存储优化
J矩阵的对称性和稀疏性带来优化空间：
对称性压缩
只存储上三角元素，节省近50%存储空间。计算时：if (i > j) J_val = J_mem[j][i]; // 读取转置位置 else J_val = J_mem[i][j];块状存储
将矩阵划分为16x16子块，每个块连续存储。这提高突发传输效率，降低存储器访问延迟。稀疏矩阵处理
对零元素占比高的场景，采用COO格式存储非零元素：typedef struct { uint8_t row; uint8_t col; int8_t val; } J_entry;
4.2 加热机制实现
加热项$\gamma y_i(t_n)\Delta t$的硬件实现需注意：
随机数生成
采用32位LFSR产生伪随机序列，经缩放后得到$\gamma$值：always @(posedge clk) begin lfsr <= {lfsr[30:0], lfsr[31] ^ lfsr[21] ^ lfsr[1]}; gamma = $signed(lfsr[15:0]) * GAMMA_SCALE; end动量更新
TE单元中的更新逻辑修改为：y_new = y_old + ( (a0 - a)*x + c0*acc )*dt + gamma*y_old*dt;
4.3 参数化设计
通过SystemVerilog参数实现灵活配置：
module dSB_core #( parameter N_SPIN = 256, parameter PC = 16, parameter PR = 16, parameter PB = 4, parameter FIXED_WIDTH = 10, parameter FIXED_FRAC = 7 ) ( input clk, input rst_n, // ...其他接口 );这种设计允许用户根据目标FPGA资源调整：
- 并行度参数（Pc/Pr/Pb）
- 问题规模（N_SPIN）
- 数值精度（FIXED_WIDTH/FIXED_FRAC）
5. 实现结果与验证
5.1 资源利用率
在KV260（Artix-7架构）上的资源占用：
| 资源类型 | 使用量 | 总量 | 利用率 |
|---|---|---|---|
| LUT | 28,421 | 154,000 | 18% |
| FF | 36,752 | 308,000 | 12% |
| DSP48 | 48 | 360 | 13% |
| BRAM | 12 | 416 | 3% |
配置参数：Pc=16, Pr=16, Pb=4, 10-bit定点
5.2 性能指标
求解256-spin max-cut问题（G-set G7）的性能：
| 指标 | 数值 |
|---|---|
| 时钟频率 | 250 MHz |
| 单次迭代周期数 | 68 cycles |
| 总迭代步数 | 10,000 |
| 总计算时间 | 2.72 ms |
| 能量效率 | 12 pJ/spin/step |
相比同平台上的bSB实现，dSB展现出：
- 速度提升1.8倍
- 精度提高15%
- 功耗降低22%
5.3 验证案例
我们测试了两个经典问题：
Max-cut问题
使用G-set标准测试集，在G7图（800节点）上的切割值达到2000，接近已知最优解（2016）。背包问题
100物品实例中，最优解发现率78%，平均误差<2%。图4展示了不同算法随迭代步数的收敛曲线。
6. 使用指南与扩展建议
6.1 快速部署步骤
环境准备
git clone https://github.com/your_repo/dSB-FPGA cd dSB-FPGA vivado -mode batch -source scripts/synth_kv260.tcl参数配置
修改include/params.svh中的关键参数：`define N_SPIN 256 `define PC 16 `define USE_HEATING 1问题输入
准备J矩阵和h向量（CSV格式），通过Python脚本转换：from utils import convert_problem convert_problem("maxcut_G7.csv", format="dSB")
6.2 常见问题排查
收敛速度慢
- 检查$c_0$是否按$\lambda_{MAX}$正确设置
- 尝试增大$\gamma$值（0.3~0.7范围）
- 确保初始$x_i$在$[-0.1,0.1]$随机分布
结果不理想
- 增加$n_{steps}$到20000以上
- 尝试多次运行取最优解
- 验证J矩阵和h向量的符号是否正确
时序违例
- 降低时钟频率10%~20%
- 减少Pc/Pr参数值
- 插入流水线寄存器
6.3 扩展方向
混合精度计算
对J矩阵采用8-bit，内部累加用16-bit，平衡精度和资源。动态加热策略
随迭代步数衰减$\gamma$值，模拟退火效果： $$ \gamma(t) = \gamma_0 \cdot e^{-t/\tau} $$多FPGA扩展
通过高速串行链路（如Aurora）连接多块FPGA，使用图分割算法分配计算任务。
这个开源项目为组合优化提供了灵活的硬件加速方案。通过调整并行度参数，可以在各类FPGA平台上实现从边缘计算到数据中心的部署。我们期待社区共同完善这一架构，推动伊辛机器在实际应用中的发展。