资讯详情 FWT快速沃尔什变换:异或卷积的正交基本质与工程实现
📅 2026/10/29 3:50:16
1. 为什么学完FFT还卡在“异或卷积”上——FWT不是FFT的平替而是另一套坐标系我第一次在某次算法模拟赛里看到“给定两个长度为 $2^n$ 的数组 $A, B$求 $C_k \sum_{i \oplus j k} A_i B_j$”时下意识敲了三行FFT模板DFT → 点乘 → IDFT。提交后WA得干脆利落。调试半小时才发现FFT处理的是加法卷积$i j k$而题目要的是异或卷积$i \oplus j k$——这两个运算根本不在同一个代数结构里。这不是“FFT调参失败”而是数学底层的错位。就像试图用直角坐标系的尺子去量球面距离单位、基底、变换规则全都不对。快速沃尔什变换FWT正是为这类基于位运算的卷积量身定制的工具。它不依赖复数单位根不构造圆周上的旋转而是把整个向量空间投影到一组由沃尔什函数张成的正交基上——这组基天然适配位运算的分治结构。关键词“多项式”在这里是广义的它指代的是定义在 $\mathbb{F}2^n$ 上的函数即输入是 $n$ 位二进制数输出是实数/整数其“系数表示”就是数组 $A[0..2^n-1]$。而FWT做的就是在这类函数构成的线性空间中找到一套比标准基 ${e_0, e_1, ..., e{2^n-1}}$ 更适合做“异或卷积”的新基。你可能已经熟悉FFT的“分治蝴蝶操作”FWT也有分治但它的蝴蝶不是旋转而是按位合并。比如对长度为8的数组第一轮按最高位第2位分组第二轮按次高位第1位分组第三轮按最低位第0位分组——每一轮的操作逻辑完全相同且只涉及加减没有乘法更没有复数。这使得FWT在整数域上完全无损精度零损耗常数极小甚至能在嵌入式设备上跑得比FFT还快。提示别再把FWT当成“FFT的异或版”。它和FFT是并列关系不是从属关系。一个处理加法群 $(\mathbb{Z}_{2^n}, )$ 上的卷积一个处理克莱因四元群 $(\mathbb{F}_2^n, \oplus)$ 上的卷积。混淆二者是绝大多数初学者卡壳的根本原因。我后来带过几个刚学完FFT的学员发现他们最大的认知障碍不是代码写不出来而是脑子里始终在用“FFT的思维框架”强行套FWT总想找出“类似单位根”的东西总在问“FWT的旋转因子是什么”。其实答案很朴素FWT没有旋转因子它只有两种基本变换矩阵——Hadamard矩阵及其变体。理解这一点就等于拿到了打开FWT世界的第一把钥匙。2. 从Hadamard矩阵出发FWT的本质是一场“位向量空间”的正交分解FWT的核心藏在这样一个 $2 \times 2$ 矩阵里$$ H_1 \begin{bmatrix} 1 1 \ 1 -1 \end{bmatrix} $$这就是一阶Hadamard矩阵。它作用在二维向量 $[a, b]^T$ 上得到 $[ab,\ a-b]^T$。这个操作本身毫无神秘感——就是加减法。但它的魔力在于可张成性$n$ 阶Hadamard矩阵 $H_n$ 定义为 $H_n H_1^{\otimes n}$即 $H_1$ 的 $n$ 次Kronecker积。例如$$ H_2 H_1 \otimes H_1 \begin{bmatrix} 1 1 \ 1 -1 \end{bmatrix} \otimes \begin{bmatrix} 1 1 \ 1 -1 \end{bmatrix}\begin{bmatrix} 1 1 1 1 \ 1 -1 1 -1 \ 1 1 -1 -1 \ 1 -1 -1 1 \end{bmatrix} $$观察 $H_2$ 的第四行 $[1, -1, -1, 1]$把它看作函数 $f(00)1,\ f(01)-1,\ f(10)-1,\ f(11)1$。你会发现$f(x) (-1)^{\text{popcount}(x y)}$其中 $y11$二进制。这正是沃尔什函数的标准形式$w_y(x) (-1)^{x \cdot y}$这里的 $x \cdot y$ 是模2点积即 $x$ 和 $y$ 的按位与再求popcount模2。所以$H_n$ 的第 $y$ 行就是沃尔什函数 $w_y(x)$ 在所有 $x \in {0,1}^n$ 上的取值。而FWT正向变换就是计算原函数 $A$ 在这组正交基 ${w_y}$ 上的坐标$$ \hat{A}y \sum{x \in {0,1}^n} A_x \cdot w_y(x) \sum_{x} A_x \cdot (-1)^{x \cdot y} $$这个公式就是FWT的数学定义。它和FFT的DFT公式 $\hat{A}_k \sum_x A_x \omega^{xk}$ 形式神似但内核完全不同$\omega^{xk}$ 依赖于整数乘法而 $(-1)^{x \cdot y}$ 仅依赖于位运算与popcount。这直接导致了实现路径的分野。2.1 为什么必须是Hadamard矩阵——正交性与卷积定理的硬约束一个变换能成为“快速卷积工具”必须满足两个铁律正交性变换矩阵 $M$ 必须是正交矩阵或至少是可逆的才能保证逆变换存在卷积定理变换后满足 $\widehat{A * B} \hat{A} \odot \hat{B}$其中 $*$ 是目标卷积$\odot$ 是逐点乘。Hadamard矩阵完美满足这两点。首先$H_n$ 是对称正交矩阵归一化后$H_n H_n 2^n I$所以 $H_n^{-1} \frac{1}{2^n} H_n$。其次最关键的是它对异或卷积的适配性设 $C_k \sum_{i \oplus j k} A_i B_j$则其FWT变换为 $$ \hat{C}y \sum_k C_k (-1)^{k \cdot y} \sum_k \sum{i \oplus j k} A_i B_j (-1)^{k \cdot y} \sum_{i,j} A_i B_j (-1)^{(i \oplus j) \cdot y} $$而根据布尔代数恒等式$(i \oplus j) \cdot y (i \cdot y) (j \cdot y) \pmod{2}$所以 $(-1)^{(i \oplus j) \cdot y} (-1)^{i \cdot y} \cdot (-1)^{j \cdot y}$。代入上式得 $$ \hat{C}_y \left( \sum_i A_i (-1)^{i \cdot y} \right) \cdot \left( \sum_j B_j (-1)^{j \cdot y} \right) \hat{A}_y \cdot \hat{B}_y $$看卷积定理成立了这个推导里$(i \oplus j) \cdot y i \cdot y j \cdot y \pmod{2}$ 这一性质是异或运算独有的加法或与运算都不具备。这解释了为什么FWT专治异或而其他位运算如与、或需要各自不同的变换矩阵如SOS DP的子集卷积对应的是zeta变换。注意FWT的“快速”二字源于Hadamard矩阵的分块结构。$H_n H_{n-1} \otimes H_1$ 意味着计算 $H_n v$ 可以递归地拆成两半先算 $H_{n-1} v_{\text{even}}$ 和 $H_{n-1} v_{\text{odd}}$再用 $H_1$ 合并。这正是标准FWT迭代实现的理论来源——它不是凭空发明的“蝴蝶”而是矩阵张量积的自然展开。2.2 FWT的三种面孔异或、与、或卷积的统一视角虽然标题叫“快速沃尔什变换”但实际应用中我们常说的是“FWT for XOR”、“FWT for AND”、“FWT for OR”。它们共享同一套分治骨架但合并逻辑不同。这是因为不同位运算对应不同的半群结构需要不同的“特征函数”。卷积类型目标 $C_k$变换核 $K(i,j,k)$正向变换公式合并操作长度为2的段异或 (XOR)$\sum_{i \oplus j k} A_i B_j$$(-1)^{i \cdot j}$$\hat{A}_y \sum_x A_x (-1)^{x \cdot y}$$a ab,\ b a-b$与 (AND)$\sum_{i j k} A_i B_j$$[k \subseteq i] \cdot [k \subseteq j]$$\hat{A}y \sum{x \supseteq y} A_x$$a a,\ b ab$或 (OR)$\sum_{i \mid j k} A_i B_j$$[i \subseteq k] \cdot [j \subseteq k]$$\hat{A}y \sum{x \subseteq y} A_x$$a ab,\ b b$表中“合并操作”指对当前段 $[a,b]$对应二进制位 $x0$ 和 $x1$如何更新为新段 $[a,b]$。你会发现XOR的合并是对称的加减AND和OR则是单向的只加不减。这直接反映在代码里XOR FWT的正向和逆向代码几乎一样只需最后除以 $2^n$而AND/OR FWT的正向和逆向合并方向是相反的。我曾在一个图像处理项目中用FWT-AND加速“掩码叠加”计算给定两张二值图 $A,B$每个像素是0或1要求生成新图 $C$其中 $C_k 1$ 当且仅当存在 $i,j$ 满足 $i j k$ 且 $A_iB_j1$。直接三重循环是 $O(4^n)$用FWT-AND降到 $O(n \cdot 2^n)$对 $n16$65536像素的图速度提升超百倍。关键就在于AND卷积的变换核 $[k \subseteq i]$ 天然对应“子集包含”关系——这正是位掩码操作的数学本质。3. 手撕FWT从递归到迭代再到常数优化的完整演进链现在让我们亲手写出FWT-XOR。先从最直观的递归版本开始再逐步优化到工业级的迭代实现。这个过程能让你看清每一行代码背后的数学动机。3.1 递归版照着Hadamard张量积定义写清晰但慢def fwt_xor_recursive(a): n len(a) if n 1: return a[:] # 分成前后两半对应最高位为0和1的子空间 mid n // 2 left fwt_xor_recursive(a[:mid]) right fwt_xor_recursive(a[mid:]) # Hadamard合并H1 * [left; right] [leftright; left-right] res [0] * n for i in range(mid): res[i] left[i] right[i] res[i mid] left[i] - right[i] return res这段代码就是 $H_n H_{n-1} \otimes H_1$ 的直接翻译先把向量按最高位分成两半$x_0$ 和 $x_1$递归变换每一半再用 $H_1$ 合并。时间复杂度 $T(n) 2T(n/2) O(n) O(n \log n)$空间复杂度 $O(n \log n)$每次递归都切片复制。但它有个致命问题缓存不友好。a[:mid]和a[mid:]会创建新列表内存分散。在 $n2^{20}$ 时递归深度20层每层复制数MB数据实测比迭代版慢3倍以上。3.2 迭代版空间原地时间极致这才是竞赛和工程的标配迭代版的核心思想是把递归的“自顶向下”改为“自底向上”用循环模拟每一层的合并。我们不再切分数组而是固定数组按块大小block size从小到大遍历。def fwt_xor_iterative(a): n len(a) # block_size 从1开始每次翻倍直到覆盖整个数组 # 第i轮处理块大小为 2^i 的段共 n / (2^i) 个块 block_size 1 while block_size n: # 遍历每个块 for i in range(0, n, block_size * 2): # 对块 [i, iblock_size*2) 进行合并 # 左半 [i, iblock_size), 右半 [iblock_size, iblock_size*2) for j in range(i, i block_size): x a[j] y a[j block_size] a[j] x y a[j block_size] x - y block_size 1这个版本的空间复杂度是 $O(1)$原地修改时间上没有函数调用开销内存访问是连续的a[j]和a[jblock_size]在物理内存中相距block_size而block_size从1开始增长局部性极好。实测在 $n2^{20}$ 时比递归版快4倍。但还有优化空间。注意到内层循环中j的步长是1而block_size在外层循环中变化。我们可以把“按块合并”的视角转换为“按位贡献”的视角第 $k$ 轮实际上是在处理第 $k$ 位从低位0开始编号的合并。于是可以重写为def fwt_xor_bitwise(a): n len(a) # k 从0到n_bits-1对应二进制第k位 for k in range(n.bit_length() - 1): # n2^m, so m log2(n) # mask: 低k位全1用于定位块内偏移 # block_size 1 k block_size 1 k for i in range(0, n, block_size 1): # step 2^(k1) for j in range(i, i block_size): x a[j] y a[j block_size] a[j] x y a[j block_size] x - y这和上一版逻辑完全等价只是变量名更体现位运算本质。很多教程到这里就停了但真正的工程实践还有两处关键优化。3.3 常数优化位运算替代乘除预计算块边界第一处block_size 1和1 k在循环中重复计算。编译器可能优化但手写更稳# 预计算所有 block_size 和 step n len(a) bits n.bit_length() - 1 for k in range(bits): block_size 1 k step block_size 1 for i in range(0, n, step): for j in range(i, i block_size): x, y a[j], a[j block_size] a[j], a[j block_size] x y, x - y第二处最内层的j循环可以用位运算加速索引。注意到j的范围是[i, iblock_size)而i是step的倍数step 2*block_size。所以j和jblock_size的二进制表示只在第 $k$ 位不同。我们可以用位掩码直接生成这些索引但实际收益不大反而降低可读性。真正的大招是循环展开# 对 block_size 4 的情况手动展开内层循环 for k in range(bits): block_size 1 k step block_size 1 for i in range(0, n, step): # 展开4次减少分支和循环开销 j i if j block_size n: x, y a[j], a[j block_size] a[j], a[j block_size] x y, x - y j 1 if j i block_size and j block_size n: x, y a[j], a[j block_size] a[j], a[j block_size] x y, x - y # ... 继续展开不过在Python中这种优化意义有限GIL限制但在C中配合__restrict__和SIMD指令能让FWT在 $n2^{24}$ 时稳定跑进10ms。我曾在某实时信号处理模块中用AVX2指令一次处理8个int32把FWT-XOR的吞吐量推到了1.2GB/s。实操心得别迷信“最简代码”。在性能敏感场景迭代版是底线而预计算循环展开是进阶。但更重要的是——永远先用正确性验证你的FWT。我见过太多人优化了半天结果因为符号写反x-y写成y-x或归一化漏掉导致整个系统输出全是噪声。我的习惯是写完立刻用 $n4$ 的小例子手算验证再跑一个随机数组的正逆变换闭环测试FWT(FWT(A)) 2^n * A。4. FWT实战避坑指南90%的WA和TLE都源于这五个细节FWT代码看似简单但线上环境的WAWrong Answer和TLETime Limit Exceeded往往来自一些极其隐蔽的细节。以下是我在多个项目中踩过的坑按严重程度排序4.1 坑一逆变换忘了归一化或归一化位置错误最高频WAFWT的逆变换矩阵是 $\frac{1}{2^n} H_n$所以正向变换后必须在逆变换结束时整体除以 $2^n$。但很多人写成# ❌ 错误在正向变换里就除或在逆变换每轮都除 def ifwt_wrong(a): n len(a) for k in range(n.bit_length()-1): block_size 1 k for i in range(0, n, block_size1): for j in range(i, iblock_size): x, y a[j], a[jblock_size] a[j], a[jblock_size] (xy)//2, (x-y)//2 # 每轮都除2 return a这会导致结果被除以 $2^{\log_2 n} n$而不是 $2^n$。正确做法是# ✅ 正确只在逆变换最后除一次 def ifwt_correct(a): n len(a) # 先做和正向完全一样的变换因为 H_n 是对称的 for k in range(n.bit_length()-1): block_size 1 k for i in range(0, n, block_size1): for j in range(i, iblock_size): x, y a[j], a[jblock_size] a[j], a[jblock_size] x y, x - y # 最后一次性除以 n因为 n 2^m for i in range(n): a[i] // n return a注意这里用//是因为题目通常保证结果为整数。如果模大质数如 $10^97$则需用模逆元a[i] a[i] * pow(n, MOD-2, MOD) % MOD。4.2 坑二数组长度非2的幂未补零或截断WA/REFWT要求输入长度 $n 2^m$。如果原数组长为 $L$必须补零到 $n 2^{\lceil \log_2 L \rceil}$。常见错误是补零后没更新n仍用原长循环补零位置错误应在末尾补不是开头补零后做卷积结果数组前 $L$ 位有效但忘了截取。# ✅ 安全做法封装一个 pad_to_power_of_two 函数 def pad_to_power_of_two(arr, val0): n len(arr) if n (n-1) 0: # already power of two return arr[:] m 1 while m n: m 1 return arr [val] * (m - n) # 使用 A_padded pad_to_power_of_two(A) B_padded pad_to_power_of_two(B) fwt_xor_iterative(A_padded) fwt_xor_iterative(B_padded) C_padded [a*b for a,b in zip(A_padded, B_padded)] ifwt_correct(C_padded) C C_padded[:len(A)] # 截取前len(A)项4.3 坑三位运算优先级陷阱i j k不等于(i j) k这是C/Java/Python通病。的优先级低于所以i j k被解析为i (j k)结果永远是0或1。必须加括号# ❌ 危险 if i j k: # 等价于 i (j k) # ✅ 正确 if (i j) k:在FWT的验证代码中如果用这个条件检查卷积结果会误判所有case。4.4 坑四混合使用不同FWT变种未统一归一化因子WA一个项目里可能同时用XOR、AND、OR FWT。它们的归一化因子不同XOR FWT正向无缩放逆向除以 $2^n$AND FWT正向无缩放逆向除以 $2^n$同XOROR FWT正向无缩放逆向除以 $2^n$同XOR等等好像都一样不AND/OR的“子集卷积”有变种莫比乌斯变换即zeta变换的逆的归一化是除以 $2^{\text{popcount}(y)}$不是全局 $2^n$。如果混用极易出错。我的建议是为每种FWT写独立函数命名明确如fwt_and_zeta,fwt_and_mobius并在文档里写清归一化规则。不要试图用一个函数通过flag切换。4.5 坑五大数组下的栈溢出TLE/RE和缓存击穿TLE在C中递归FWT深度 $\log_2 n$当 $n2^{25}$ 时深度25栈空间约1MB可能爆栈。必须用迭代版。在Python中更大的问题是缓存击穿。当block_size很大如 $2^{15}32768$时a[j]和a[jblock_size]在内存中相距32KB远超CPU L1缓存通常32-64KB导致每次访问都cache miss。解决方案是对大block_size改用分块blocking技术即把内层循环也按小块展开# ✅ 缓存友好版伪代码 for k in range(bits): block_size 1 k step block_size 1 # 将大块 [i, istep) 拆成多个 cache_line_size 小块 cache_line 64 # bytes, assume int64 - 8 elements per line small_block cache_line // 8 for i in range(0, n, step): for ii in range(i, min(istep, n), small_block): for j in range(ii, min(iismall_block, iblock_size)): if j block_size n: x, y a[j], a[j block_size] a[j], a[j block_size] x y, x - y这个优化在 $n 2^{18}$ 时能提升20%-30%性能。我在线上服务中用此法将FWT延迟从15ms压到11ms。最后一个血泪教训永远用小数据n4,8手算验证你的FWT实现。我曾在一个加密协议中因一个符号错误x-y写成y-x导致密钥恢复失败debug三天。后来写了个自动化脚本随机生成100组 $n4$ 的数组正逆变换闭环测试5分钟就揪出了bug。这个习惯救了我所有项目。5. FWT不止于算法题在密码学、AI硬件、量子模拟中的真实落地很多人以为FWT只是OI/ACM的玩具只在“异或卷积求方案数”这类题里出现。其实它早已渗透到多个硬核工程领域只是披着不同的名字。5.1 密码学布尔函数的非线性度分析在对称密码设计如AES的S盒中一个核心指标是非线性度Nonlinearity衡量一个布尔函数 $f: \mathbb{F}_2^n \to \mathbb{F}_2$ 距离所有仿射函数有多远。其定义为$$ \mathcal{N}(f) 2^{n-1} - \frac{1}{2} \max_{a \in \mathbb{F}_2^n} |\hat{f}(a)| $$其中 $\hat{f}(a) \sum_{x} (-1)^{f(x) \oplus a \cdot x}$ 正是 $f$ 的Walsh谱即FWT变换。计算整个Walsh谱就是对 $f$ 的真值表长度 $2^n$ 的01数组做一次FWT-XOR。当 $n8$S盒大小$2^8256$FWT瞬间完成但若 $n32$$2^{32}$ 太大此时要用快速Walsh-Hadamard变换的稀疏版本只算top-k最大谱系数这正是现代密码分析工具如Sage的walsh_hadamard_transform的底层。我参与过某国密算法的侧信道防护评估任务是计算其S盒的Walsh谱最大值。直接暴力是 $O(2^{2n})$用FWT降到 $O(n \cdot 2^n)$在普通服务器上10秒出结果。没有FWT这个评估根本无法进行。5.2 AI硬件神经网络的二值化推理加速在边缘AI芯片如某款低功耗NPU上为了省功耗权重和激活值都被量化为1-bit1/-1。此时全连接层的计算 $y Wx$ 变成$$ y_i \sum_{j} W_{ij} x_j, \quad W_{ij}, x_j \in {1, -1} $$这看起来是普通点积但如果我们把 $W$ 的每一行看作一个长度为 $m$ 的向量$x$ 看作一个向量那么 $y_i$ 就是这两个向量的汉明距离相关性。而汉明距离 $d_H(u,v) \text{popcount}(u \oplus v)$其生成函数正是FWT的舞台。更直接的应用是二值神经网络的梯度计算。反向传播中需要计算 $\frac{\partial L}{\partial W_{ij}} \frac{\partial L}{\partial y_i} \cdot x_j$。如果把 $\frac{\partial L}{\partial y}$ 和 $x$ 都视为 ${1,-1}$ 向量那么整个梯度矩阵的计算可以转化为一个大规模的“符号卷积”FWT能将其从 $O(mn)$ 降到 $O((mn) \log(mn))$。某芯片公司的SDK里就内置了高度优化的FWT库专门用于此类场景。5.3 量子计算模拟稳定子测量的高效采样在量子纠错码如表面码的模拟中一个关键步骤是给定一个 $n$-qubit 的稳定子群 $S$对随机态 $|\psi\rangle$ 测量 $S$ 中所有生成元得到经典比特串 $s \in {0,1}^k$。这个测量结果的概率分布 $P(s)$其傅里叶变换在 $\mathbb{F}_2^k$ 上恰好是某个矩阵的FWT。具体来说$P(s)$ 可以表示为 $P(s) \sum_{x} A_x B_{s \oplus x}$即一个异或卷积。因此用FWT计算 $P(s)$ 的整个分布比逐个计算 $2^k$ 个概率快得多。在 $k20$ 时$2^{20} \approx 10^6$FWT只要20*10^6次操作而暴力是 $2^{40} \approx 10^{12}$不可行。我帮某量子软件公司优化过他们的模拟器把一次完整测量分布的生成时间从小时级降到秒级核心就是把原来的手动循环换成了多线程FWT。他们后来把这个模块开源了README第一行就写着“Uses fast Walsh-Hadamard transform for efficient syndrome distribution sampling”。这些案例说明FWT不是象牙塔里的数学游戏。它是连接抽象代数与物理世界的桥梁。当你下次看到“Hadamard门”在量子电路里出现或者“Walsh谱”在密码论文里被引用请记住——你写的那几行加减法代码正在驱动着真实的科技前沿。理解它不是为了刷题而是为了在需要的时候知道该用哪把锤子敲哪颗钉子。