HomeArchiveBlog


Original contents are licensed under CC BY-NC 4.0. All rights reserved © 2026 Kai.
Back to Archives
Advanced FPGA HLS: Non-linear Programming

Explore advanced non-linear programming techniques for FPGA High-Level Synthesis (HLS) to optimize performance and resource utilization.

Tue Jan 13 2026
Tue Jan 13 2026
FPGAHLSNon-Linear Programming
On this page
  • Advanced FPGA HLS: Non-Linear Programming
    • 算法示例
    • 算法分析
      • 循环展开
      • 循环流水化
      • 循环分块
      • 循环交换
      • 资源约束
      • 完整模型
    • 考虑 BRAM 约束
    • 考虑全局内存访问
    • 考虑数据流优化
    • 参考文献

Advanced FPGA HLS: Non-Linear Programming

FPGA HLS 的优化过程实际上是一个组合优化问题, 涉及多种资源约束, 例如 DSPs, BRAMs 约束; 以及多种性能目标, 例如全局总延迟, 吞吐量等. 对于组合优化问题, 数学上最直接的方法是使用线性规划(LP) 或者非线性规划 (NLP) 直接建模求解, 并且理论上能找到全局最优解 (尽管大多数情况下不能在合理时间内找到). 在调度过程中 LP 在工业界广为使用, 它适用于细粒度的调度问题 (将每条操作实例映射到一个时间戳上), 适合已经给定一段操作序列, 然后求解最优解. 如果从更高的算法层次来考虑问题, 例如循环分块策略, 循环展开因子等, 则问题通常不能简单用 LP 建模, 转而需要使用 NLP 方法. 本文介绍 FPGA HLS 中如何使用 NLP 方法求解最优循环优化策略.

算法示例

我们实现一个简单的 3MM (Triple Matrix Multiplication) 算法, 该算法包含三个矩阵乘法操作, 数学表达为:

E=AB,F=CD,G=EFE = AB, \quad F = CD, \quad G = EFE=AB,F=CD,G=EF

代码实现如下

template<typename T>
void triple_mm(T A[16][16], T B[16][16], T C[16][16], T D[16][16], T E[16][16], T F[16][16], T G[16][16]) {
    // First Matrix Multiplication E = A * B
    for (int i = 0; i < 16; i++) {
        for (int j = 0; j < 16; j++) {
            E[i][j] = 0;
            for (int k = 0; k < 16; k++) {
                E[i][j] += A[i][k] * B[k][j];
            }
        }
    }
    // Second Matrix Multiplication F = C * D
    for (int i = 0; i < 16; i++) {
        for (int j = 0; j < 16; j++) {
            F[i][j] = 0;
            for (int k = 0; k < 16; k++) {
                F[i][j] += C[i][k] * D[k][j];
            }
        }
    }
    // Third Matrix Multiplication G = E * F
    for (int i = 0; i < 16; i++) {
        for (int j = 0; j < 16; j++) {
            G[i][j] = 0;
            for (int k = 0; k < 16; k++) {
                G[i][j] += E[i][k] * F[k][j];
            }
        }
    }
}

这里面有一个简单的数据流依赖, 即第三个矩阵乘法 G=EFG = EFG=EF 依赖于前两个矩阵乘法的结果 EEE 和 FFF. 我们希望通过调整每个矩阵乘法的循环优化策略 (例如循环展开因子, 循环分块大小等) 来最小化整体延迟, 同时满足 FPGA 资源约束 (例如 DSP 数量).

算法分析

循环展开

在最简单的情况下, 循环没有流水线, 没有展开, 那么每个矩阵乘法的延迟即为

L=16×16×16×(tmul+tadd)L = 16 \times 16 \times 16 \times (t_\text{mul} + t_\text{add})L=16×16×16×(tmul​+tadd​)

此时我们假设加载数据和写回数据的延迟可以忽略不计. 整个模块的总延迟即为三个矩阵乘法延迟之和 3L3L3L.

由于矩阵乘法的 k 维度上存在数据依赖, 如果对 k 维度进行展开会触发树形归约优化 (Tree Reduction), 使得 k 上的操作可以部分并行执行, 如下图所示

假设 k 维度的展开因子为 uku_kuk​, 且 uku_kuk​ 能整除 16, 那执行一次完整的 k 循环需要的周期数为

Lk=16uk×(tmul+⌈log⁡2uk⌉×tadd)L_k = \frac{16}{u_k} \times (t_\text{mul} + \lceil \log_2{u_k} \rceil \times t_\text{add})Lk​=uk​16​×(tmul​+⌈log2​uk​⌉×tadd​)

总周期数为

L=16×16×LkL = 16 \times 16 \times L_kL=16×16×Lk​

如果不在 k 上做展开, 由于 i, j 和 k 之间的顺序可以交换, 我们可以将 k 循环挪到外面, 将 i 或者 j 挪到内层, 将 i 和 j 进行展开. 由于 i 和 j 维度上都没有数据依赖, 因此这两个维度上展开后的总延迟即为

Li=Lj=16u×(tmul+tadd)L_i = L_j = \frac{16}{u} \times (t_\text{mul} + t_\text{add})Li​=Lj​=u16​×(tmul​+tadd​)

总延迟为

L=16×16×Li=16×16×LjL = 16 \times 16 \times L_i = 16 \times 16 \times L_jL=16×16×Li​=16×16×Lj​

因此循环展开对总延迟的影响, 如果在展开维度上存在数据依赖并启用树形归约, 则展开因子对延迟的影响会引入一个非线性项 ⌈log⁡2uk⌉\lceil \log_2{u_k} \rceil⌈log2​uk​⌉.

循环流水化

现在假设在 k 维度上展开因子为 uku_kuk​, 在 i, j 维度上均不做展开处理. 我们选择 i 和 j 中的其中一个维度进行流水化. 由于在 i 和 j 维度上并不存在数据依赖, 因此在这两个维度上的启动间隔理论上可以达到 1 (暂时不考虑资源约束). 那么在流水化维度上执行一次完整的循环需要的周期数为

Li=Lj=(16−1)×1+LkL_i = L_j = (16 - 1) \times 1 + L_kLi​=Lj​=(16−1)×1+Lk​

总延迟为

L=16×Li=16×LjL = 16 \times L_i = 16 \times L_jL=16×Li​=16×Lj​

由于 i 和 j 只有一个维度可以被流水化 (通常在 HLS 中只能对一个循环进行流水化), 因此引入决策变量 pi,pjp_i, p_jpi​,pj​, 它们都是二元变量, 且满足约束条件

pi+pj=1p_i + p_j = 1pi​+pj​=1

表示必须且只能选择其中一个维度进行流水化. 那么总延迟可以表示为

L=16×(pi×Li+pj×Lj)L = 16 \times (p_i \times L_i + p_j \times L_j)L=16×(pi​×Li​+pj​×Lj​)

至此我们就利用循环展开因子 uku_kuk​ 和流水化决策变量 pi,pjp_i, p_jpi​,pj​ 将这两种优化方式纳入到延迟模型中.

循环分块

接着考虑如何将循环分块纳入到模型中. 假设在每个维度上都分成点循环和块循环两部分, 例如在 i 维度上分块大小为 i0i_0i0​, i1i_1i1​, 那么有约束条件

i0×i1=16,j0×j1=16,k0×k1=16i_0 \times i_1 = 16, j_0 \times j_1 = 16, k_0 \times k_1 = 16i0​×i1​=16,j0​×j1​=16,k0​×k1​=16

我们采用一种微架构设计模式, 将内层的点循环视作 PE Cluster, 由一群相同且并行的基本功能单元组成; 外层的块循环视作流水化的任务流 (Task Pipeline). 具体实现上也就是内层的所有点循环都会被完全展开, 形成一个个独立的 PE Cluster, 然后外层块循环中的某个维度会被流水化, 形成一个任务流水线. 例如

for (int i0 = 0; i0 < 4; i0++) {
    for (int j0 = 0; j0 < 4; j0++) { // task pipeline loop
        #pragma HLS pipeline ii=1
        for (int k0 = 0; k0 < 4; k0++) {
            // point loop unroll
            for (int i1 = 0; i1 < 4; i1++) {
                #pragma HLS unroll
                for (int j1 = 0; j1 < 4; j1++) {
                    #pragma HLS unroll
                    for (int k1 = 0; k1 < 4; k1++) {
                        #pragma HLS unroll
                        // compute
                    }
                }
            }
        }
    }
}

具体的并行程度和流水线深度通过循环分块大小来决定. 内层所有点循环的延迟通过前面的展开模型来计算, 由于完全展开, 因此对于有或没有数据依赖的维度, 延迟分别为

Li1=Lj1=tmul+taddL_{i_1} = L_{j_1} = t_\text{mul} + t_\text{add}Li1​​=Lj1​​=tmul​+tadd​ Lk1=tmul+⌈log⁡2u⌉×taddL_{k_1} = t_\text{mul} + \lceil \log_2{u} \rceil \times t_\text{add}Lk1​​=tmul​+⌈log2​u⌉×tadd​

点循环的总延迟应该取决于最大的那个内层循环延迟, 即

Lpoint=max⁡(Li1,Lj1,Lk1)L_\text{point} = \max(L_{i_1}, L_{j_1}, L_{k_1})Lpoint​=max(Li1​​,Lj1​​,Lk1​​)

外层块循环的延迟则通过流水线模型进行计算, 对于被流水化的维度, 延迟为

Lp=(N−1)×II+LinnerL_p = (N - 1) \times II + L_\text{inner}Lp​=(N−1)×II+Linner​

其中 NNN 为当前维度的块循环大小, II 为启动间隔, LinnerL_\text{inner}Linner​ 表示所有在它内层的循环延迟, 而不仅仅是点循环的延迟. 块循环中只能挑一个维度进行流水化, 因此我们引入决策变量 pi,pj,pkp_i, p_j, p_kpi​,pj​,pk​, 它们都是二元变量, 且满足约束条件

pi+pj+pk=1p_i + p_j + p_k = 1pi​+pj​+pk​=1

在 i 和 j 维度上, 由于没有数据依赖, 启动间隔理论上可以达到 1. 而在 k 维度上由于存在 RAW 依赖, 在树形归约下流水线稳态的启动间隔应该等于完成单次加法需要的时钟数量, 假设需要 3 个周期. 那么在当前的循环排布下, 在三个维度上分别的延迟为

Lk0=pk×((k0−1)×3+Lpoint)+(1−pk)×k0×LpointL_{k_0} = p_k \times ((k_0 - 1) \times 3 + L_\text{point}) + (1 - p_k) \times k_0 \times L_\text{point}Lk0​​=pk​×((k0​−1)×3+Lpoint​)+(1−pk​)×k0​×Lpoint​ Lj0=pj×((j0−1)×1+Lk0)+(1−pj)×j0×Lk0L_{j_0} = p_j \times ((j_0 - 1) \times 1 + L_{k_0}) + (1 - p_j) \times j_0 \times L_{k_0}Lj0​​=pj​×((j0​−1)×1+Lk0​​)+(1−pj​)×j0​×Lk0​​ Li0=pi×((i0−1)×1+Lj0)+(1−pi)×i0×Lj0L_{i_0} = p_i \times ((i_0 - 1) \times 1 + L_{j_0}) + (1 - p_i) \times i_0 \times L_{j_0}Li0​​=pi​×((i0​−1)×1+Lj0​​)+(1−pi​)×i0​×Lj0​​

那么总延迟即为最外层块循环的延迟

L=Li0L = L_{i_0}L=Li0​​

至此我们引入了循环展开因子, 循环分块大小和流水化决策变量, 并将它们纳入到延迟模型中.

循环交换

接下来我们考虑循环交换对延迟的影响. 以上循环共有 3 个维度, 并且全部是可交换的 (有些情况并不如此). 因此共有 3!=63!=63!=6 种循环排列方式. 我们引入 6 个二元决策变量 s1,s2,s3,s4,s5,s6s_1, s_2, s_3, s_4, s_5, s_6s1​,s2​,s3​,s4​,s5​,s6​, 分别对应 6 种循环排列方式, 并且满足约束条件

s1+s2+s3+s4+s5+s6=1s_1 + s_2 + s_3 + s_4 + s_5 + s_6 = 1s1​+s2​+s3​+s4​+s5​+s6​=1

表示只能选择其中一种循环排列方式. 那么最终的总延迟可以表示为

L=∑i=16si×LiL = \sum_{i=1}^{6} s_i \times L_iL=i=1∑6​si​×Li​

其中 LiL_iLi​ 表示第 iii 种循环排列方式下的总延迟. 例如在 i, j, k 顺序下的总延迟就是刚刚计算的 Li0L_{i_0}Li0​​. 不难看出循环交换实际上只影响块循环的表达式, 因为相对顺序改变会影响 LinnerL_\text{inner}Linner​ 的计算方式. 但对于点循环的延迟计算没有任何影响.

资源约束

到此我们就有了完整的, 涵盖多种循环优化方式的延迟模型. 接下来我们考虑 FPGA 资源约束. 暂时只考虑 DSP 约束. 由于资源共享的问题, 在综合之前要想获得准确估计是比较困难的. 我们这里采用悲观估计, 也就是将所有被展开的乘法操作都视作独立的 DSP 资源来计算. 那么 DSP 资源消耗可以表示为

DSP=i1×j1×(k1×3+⌈k12⌉×2)\text{DSP} = i_1 \times j_1 \times \left(k_1 \times 3 + \left\lceil \frac{k_1}{2} \right\rceil \times 2\right)DSP=i1​×j1​×(k1​×3+⌈2k1​​⌉×2)

这样考量的原因是, 一共有 i1×j1i_1 \times j_1i1​×j1​ 棵归约树, 不同的树之间资源不能共享, 并行的操作之间资源也不能共享, 但有先后顺序的操作之间资源可以共享. 每棵树有 k1k_1k1​ 个并行乘法操作, 每个乘法操作消耗 3 个 DSP (这是假设值). 由于加法树有先后执行顺序, 需要的资源取决于最宽的那一层, 计算为 ⌈k12⌉\lceil \frac{k_1}{2} \rceil⌈2k1​​⌉, 并且每个加法操作消耗 2 个 DSP (假设值). 可以观察到, 资源约束和外层块循环大小无关, 只和内层点循环相关.

最终我们希望满足资源约束

DSP≤MAX_DSP\text{DSP} \leq \text{MAX\_DSP}DSP≤MAX_DSP

实际上这个建模一定会大于实际用量, 最显然的原因就是同一归约树中的乘法没必要同时开始, 只要保证加法树的输入在需要的时候可用即可 (我们这里是基于 ASAP 调度的基本估计). 这会导致最终求解结果比较保守.

完整模型

到此我们就有了完整的 NLP 模型, 包括模型参数:

  • 基础操作延迟: tmul,taddt_\text{mul}, t_\text{add}tmul​,tadd​
  • 计算参数: 矩阵大小 16
  • 资源约束: MAX_DSP
  • 数据依赖分析结果: k 维度上的依赖链条 (由此推断出 II)

基础变量

{i0,i1,j0,j1,k0,k1,pi,pj,pk,s1,s2,s3,s4,s5,s6}\{i_0, i_1, j_0, j_1, k_0, k_1, p_i, p_j, p_k, s_1, s_2, s_3, s_4, s_5, s_6\}{i0​,i1​,j0​,j1​,k0​,k1​,pi​,pj​,pk​,s1​,s2​,s3​,s4​,s5​,s6​}

执行次数约束

i0×i1=16,j0×j1=16,k0×k1=16i_0 \times i_1 = 16, j_0 \times j_1 = 16, k_0 \times k_1 = 16i0​×i1​=16,j0​×j1​=16,k0​×k1​=16

流水化决策约束

pi+pj+pk=1p_i + p_j + p_k = 1pi​+pj​+pk​=1

循环交换决策约束

s1+s2+s3+s4+s5+s6=1s_1 + s_2 + s_3 + s_4 + s_5 + s_6 = 1s1​+s2​+s3​+s4​+s5​+s6​=1

资源约束

DSP≤MAX_DSP\text{DSP} \leq \text{MAX\_DSP}DSP≤MAX_DSP

其中 DSP 可以用基础模型变量表示出来. 以及总延迟目标函数

L=∑i=16si×LiL = \sum_{i=1}^{6} s_i \times L_iL=i=1∑6​si​×Li​

其中每个 LiL_iLi​ 也都可以用基础模型变量表示出来. 最终我们希望最小化总延迟 LLL.

交给 NLP 求解器之后, 就可以得到最优的循环优化策略 (循环展开因子, 循环分块大小, 流水化维度, 循环顺序等). 并且这个结果是在给定资源约束下的全局最优解.

我们仔细考虑一下目标函数和约束条件, 它们都包含了大量的非线性项, 例如多个模型变量的乘积, 甚至在启用归约树的情况下还包含了对数项. 这使得问题变得非常复杂, 极难在合理时间内找到最优解. 实际情况下只有在对小模块进行极致优化的时候才会采用这种方法. 并且由于 NLP 求解器 (Gurobi 等) 的许可证问题, 商业使用受限制很严重.

考虑 BRAM 约束

另一个约束是 BRAM 约束. BRAM 约束分成两种, 一种是总容量约束, 另一种是端口数约束. 总容量约束比较简单, 只需要统计所有数组的总大小即可. 实际情况下总容量限制更难触发, 更容易的是端口数约束. 由于并行需求, 要通过 #pragma HLS array_partition 等指令将数组进行分区, 以提供足够的并行访问端口.

因此 BRAM 端口数约束应该可以用内层点循环的参数大小来建模. 实际处理的时候, 首先要获取每个数组的访问模式, 也就是在每次迭代上会访问数组的哪几个点. 例如对于访问模式

f(i,j,k)↦(i,k−1),f(i,j,k)↦(i,k),f(i,j,k)↦(i,k+1)f(i,j,k) \mapsto (i,k-1), f(i,j,k) \mapsto (i,k), f(i,j,k) \mapsto (i,k+1)f(i,j,k)↦(i,k−1),f(i,j,k)↦(i,k),f(i,j,k)↦(i,k+1)

在一次迭代中会访问数组的 (i,k-1), (i,k), (i,k+1) 相邻三个点. 如果内层点循环中 i 维度上分块大小为 i1i_1i1​, k 维度上分块大小为 k1k_1k1​, 那么在一次迭代中会访问的内存点数为 i1×(k1+2)i_1 \times (k_1 + 2)i1​×(k1​+2), 要想获得完全的并行度, 则需要将数组在 i 维度上分成 i1i_1i1​ 份, 在 k 维度上分成 k1+2k_1 + 2k1​+2 份, 总共需要 i1×(k1+2)i_1 \times (k_1 + 2)i1​×(k1​+2) 个端口.

因此这种情况下 BRAM 端口数约束可以表示为 (一块 BRAM 通常有 2 个端口)

i1×(k1+2)≤2×MAX_BRAMi_1 \times (k_1 + 2) \leq 2 \times \text{MAX\_BRAM}i1​×(k1​+2)≤2×MAX_BRAM

BRAM 端口数其实不是很好建模, 主要是因为涉及到访问模式的分析. 对于多面体限制下的分析相对好做, 访问模式是迭代变量的仿射函数, 给定点循环的分块大小 (也就是迭代空间), 可以求出访问的内存点的集合, 统计这个集合的大小 (通常是关于分块大小的多项式函数) 就能得到端口数需求. 对于更复杂的访问模式, 比如间接访问, 也就是在编译期完全无法预测的, 基本上没有方法相对准确地建模. 即使是多面体限制下的访问模式, 如果出现多项式函数 (虽然几乎不会特别高次), 也会导致模型变得非常复杂, 极难求解.

考虑全局内存访问

全局内存访问也是内存来源的一个重要方面. 考虑在上面的 3MM 算法中, 将加载操作和写回操作放在哪一个迭代维度上进行是一个重要的设计决策. 例如在 i 维度上进行数据传输, 数据复用的范围会比较大, 需要更多的缓存资源 (BRAM/URAM), 但可以减少全局内存访问次数, 从而降低延迟. 相反如果在 k 维度上进行数据传输, 则数据复用的范围很小 (每个元素只用一次就扔掉), 但全局内存访问次数会增加, 导致更高的延迟.

要将全局内存访问纳入到模型中, 需要先分析每个维度上一共要访问多少个内存点. 这和统计 BRAM 端口数类似, 需要先分析访问模式, 然后结合循环分块大小来统计内存点数. 区别在于 BRAM 只统计点循环内的访问点数 (因为点循环决定并行度), 而全局内存访问需要统计块循环内的访问点数, 并且在不同的循环维度上需要传输的数据量不同.

为了考虑全局内存访问, 需要引入新的决策变量, 表示在哪个维度上对哪一块数组进行数据传输. 例如对于 3MM 算法, 我们可以引入变量 dAi,dAj,dAkdA_i, dA_j, dA_kdAi​,dAj​,dAk​, 分别表示在 i, j, k 维度上对数组 A 进行数据传输 (二元变量). 并且满足约束条件

dAi+dAj+dAk=1dA_i + dA_j + dA_k = 1dAi​+dAj​+dAk​=1

对其它的数组也有类似处理. 接着使用多面体分析方法, 结合循环分块大小来统计每个数组在每个维度上需要传输的内存点数. 最终可以得到每个数组在每个维度上的传输点数. 而数据传输需要的延迟根据内存模型建模即可, 包括位宽, 总线延迟. 例如位宽为 512 位, 每次全局访问需要 4 个周期, 那么传输 1024 个点 (假设为 32 位浮点) 需要的延迟为

Lmem=1024×32512×4L_\text{mem} = \frac{1024 \times 32}{512} \times 4Lmem​=5121024×32​×4

最终总延迟需要加上所有数组的全局内存访问延迟.

考虑数据流优化

在 FPGA HLS 中, 数据流 (Dataflow) 优化是一种重要的优化手段, 通过在模块之间建立数据流通道 (插入 FIFO), 可以实现模块间的流水线执行, 提高整体吞吐量. 但建模过程涉及到数据流理论, 请参考 Stream-HLS [1] 一文. 这里不做赘述.

参考文献

  • Stream-HLS [1]
  • NLP-DSE [2]
  • Sisyphus [2]
[1]
Basalama, S. and Cong, J. 2025. Stream-HLS: Towards Automatic Dataflow Acceleration. Proceedings of the 2025 ACM/SIGDA International Symposium on Field Programmable Gate Arrays (Monterey, CA, USA, 2025), 103–114.
[2]
Pouget, S. et al. 2025. A Unified Framework for Automated Code Transformation and Pragma Insertion. (Monterey, CA, USA, 2025), 187–198.
[3]
Pouget, S. et al. 2024. Automatic Hardware Pragma Insertion in High-Level Synthesis: A Non-Linear Programming Approach. Proceedings of the 2024 ACM/SIGDA International Symposium on Field Programmable Gate Arrays (Monterey, CA, USA, 2024), 184.