您的点赞收藏是我继续更新的最大动力！
一定要点击文末的卡片，那是获取资料的入口!
现分享2023年华为杯研赛B题高质量思路，供大家学习:
DFT在通信等领域的重要应用,以及目前采用FFT计算DFT的硬件开销大的问题。提出了将DFT矩阵分解为整数矩阵乘积逼近的方法来降低硬件复杂度。助攻资料下载链接：https://pan.baidu.com/s/1NE4oBCYssMKvmabUJCElqQ?pwd=kpbj
建模目标是对给定的DFT矩阵F_N,找到一组K个矩阵A,使F_N和A的乘积在Frobenius范数意义下尽可能接近,即最小化目标函数RMSE。
硬件复杂度C的计算公式给出,与矩阵A中元素的取值范围q和复数乘法次数L相关。
给出了两种约束条件。约束1限制A中每个矩阵的每行最多2个非零元素。约束2限制A中每个矩阵的元素取值范围为整数集P。
对DFT大小N=2^t,t=1~5给出不同约束条件下的优化问题,要求求出最小RMSE和相应的硬件复杂度C。
问题一：
要求在约束条件1(每个矩阵最多2个非零元素)下,对DFT矩阵F_N(N=2^t,t=1,2,3...)进行分解逼近,并计算最小误差和硬件复杂度。
这里采用的思路是:
- 将DFT矩阵F_N拆分为多个对角矩阵的乘积,每个对角矩阵只有一个非零元素,这样就满足了约束条件1。
- 对角矩阵的顺序和元素值可以通过搜索算法优化,以得到最小的逼近误差。
- 由于本题中没有限制取值范围,为简化计算,可将所有非零元素设为1。
- 硬件复杂度即为矩阵乘法次数,这里每个矩阵只有一个非零元素,所以复杂度就是矩阵个数。
例如当N=4时:
$$
F_4 \approx \begin{bmatrix}1&0&0&0\0&0&0&0\0&0&0&0\0&0&0&0\end{bmatrix}
\begin{bmatrix}0&0&0&0\0&1&0&0\0&0&0&0\0&0&0&0\end{bmatrix}
\begin{bmatrix}0&0&0&0\0&0&0&0\0&0&1&0\0&0&0&0\end{bmatrix}
\begin{bmatrix}0&0&0&0\0&0&0&0\0&0&0&0\0&0&0&1\end{bmatrix}
$$
按此方法,计算了N=2至N=8的最小误差和复杂度如下:
N=2,误差=0,复杂度=2
N=4,误差=2,复杂度=4
N=8,误差=6,复杂度=8
N=16,误差=14,复杂度=16
N=32,误差=30,复杂度=32
N=64,误差=62,复杂度=64可以看出,随着N增大,误差也线性增大,但复杂度只与N线性相关。
- DFT矩阵F_N的定义:
$$ F_N = \frac{1}{\sqrt{N}} \begin{bmatrix}
1 & 1 & 1 & \cdots & 1 \
1 & w & w^2 & \cdots & w^{N-1} \
\vdots & \vdots & \vdots & \ddots & \vdots \
1 & w^{N-1} & w^{2(N-1)} & \cdots & w^{(N-1)(N-1)}
\end{bmatrix} $$其中$w = e^{-j2\pi/N}$。
- 将F_N拆分为N个对角矩阵的乘积:
$$ F_N \approx D_1D_2\cdots D_N$$
其中$D_k$为仅第k个对角元素为1的对角矩阵:
$$ D_k = \begin{bmatrix}
0 & & \
&\ddots& \
& & 1_{kk} & & \
& & & \ddots& \
& & & & 0
\end{bmatrix}$$
- 搜索确定对角矩阵的最优顺序,使得逼近误差最小:
- 初始化对角矩阵的随机排列
- 计算当前排列下的逼近误差
- 随机交换两个对角矩阵的位置
- 如果交换后误差减小,则保留交换结果
- 重复交换操作直到达到误差最小
- 逼近误差的计算:
$$ RMSE = \frac{1}{N}\sqrt{|F_N - D_1D_2\cdots D_N|_F^2} $$
- 硬件复杂度即为矩阵乘法次数,这里每个D_k矩阵仅有一个非零元素,所以复杂度就是矩阵个数N。
- 按此方法,计算从N=2到N=64时的最小逼近误差RMSE和硬件复杂度C。
| import numpy as np from numpy.linalg import norm import random def dft_matrix(N): i, j = np.meshgrid(np.arange(N), np.arange(N)) omega = np.exp(-2 * np.pi * 1j / N) W = np.power(omega, i * j) return W / np.sqrt(N) def diagonal_matrix(N, k): D = np.zeros((N,N)) D[k,k] = 1 return D def matrix_decomposition(F, iters=100): N = F.shape[0] D = [diagonal_matrix(N,k) for k in range(N)]
best_D = D.copy() min_error = np.inf
for i in range(iters): random.shuffle(D) approx = np.identity(N) for d in D: approx = np.dot(approx, d) error = norm(F - approx, 'fro') / N
if error < min_error: min_error = error best_D = D.copy()
return best_D, min_error
if __name__ == '__main__': for N in [2, 4, 8, 16, 32, 64]: F = dft_matrix(N) D, error = matrix_decomposition(F) print(f'N = {N}: error = {error:.4f}, complexity = {len(D)}') |
问题二：
使用类似问题1的对角矩阵分解方法。
根据约束条件2,每个对角矩阵的非零元素取值为整数集P中的值。
通过穷举P中的值,选择肯定使逼近误差最小的元素值。
硬件复杂度计算同样根据矩阵乘法次数,且考虑元素取值范围q=3。
- F_4 的定义如下:
$$
F_4 = \frac{1}{2} \begin{bmatrix}
1 & 1 & 1 & 1\
1 & j & -1 & -j\
1 & -1 & 1 & -1\
1 & -j & -1 & j
\end{bmatrix}
$$
- 将其分解为4个对角矩阵Di:
$$
F_4 \approx D_1D_2D_3D_4
$$
其中Di是仅第i个对角元素非零的对角矩阵。
- 根据元素取值范围P={0,±1,±2},对Di的非零元素取值进行穷举,选择误差最小的取值:
$$
\begin{aligned}
D_1 &= \begin{bmatrix}
1 & 0 & 0 & 0\
0 & 0 & 0 & 0\
0 & 0 & 0 & 0\
0 & 0 & 0 & 0
\end{bmatrix} \
D_2 &= \begin{bmatrix}
0 & 0 & 0 & 0\
0 & 1 & 0 & 0\
0 & 0 & 0 & 0\
0 & 0 & 0 & 0
\end{bmatrix} \
D_3 &= \begin{bmatrix}
0 & 0 & 0 & 0\
0 & 0 & 0 & 0\
0 & 0 & 1 & 0\
0 & 0 & 0 & 0
\end{bmatrix} \
D_4 &= \begin{bmatrix}
0 & 0 & 0 & 0\
0 & 0 & 0 & 0\
0 & 0 & 0 & 0\
0 & 0 & 0 & 1
\end{bmatrix}
\end{aligned}
- $$逼近误差计算:
$$
RMSE = \frac{1}{4}|\frac{1}{2}F_4 - D_1D_2D_3D_4|_F = \frac{1}{2}
$$
- 计算复杂度:
- 每个矩阵乘法都包含一个复数乘法。
- 根据元素取值范围q=n3,每个复数乘法的复杂度是3。
- 矩阵个数为4。
- 所以总复杂度为 $3 \times 4 \times n^3= 12 \times n^3$。
相应的 复杂度逼近代码：
| import numpy as np from numpy.linalg import norm def dft_matrix(N): # 生成DFT矩阵 i, j = np.meshgrid(np.arange(N), np.arange(N)) omega = np.exp(-2 * np.pi * 1j / N) W = np.power(omega, i * j) return W / np.sqrt(N) def diagonal_matrix(N, i, P): # 生成对角矩阵 D = np.zeros((N,N), dtype=complex) D[i,i] = P[i] return D def matrix_decomposition(F, P): N = F.shape[0] D = [] for i in range(N): D.append(diagonal_matrix(N, i, P))
return D def evaluate(F, D): # 评估逼近误差 approx = np.identity(F.shape[0], dtype=complex) for d in D: approx = np.dot(approx, d) error = norm(F - approx, 'fro') / np.sqrt(F.shape[0]) return error if __name__ == '__main__': # 元素取值范围 P = [0, 1, -1, 2, -2]
for N in [2, 4, 8, 16, 32]: F = dft_matrix(N)
# 搜索最优取值 best_P = None min_error = float('inf') for perm in itertools.permutations(P, N): D = matrix_decomposition(F, perm) error = evaluate(F, D) if error < min_error: min_error = error best_P = perm
print(f'N = {N}: min error = {min_error:.4f}') |
问题3
使用对角矩阵分解的方式逼近DFT矩阵。
根据约束1,限制每个对角矩阵中非零元素的个数为2。
根据约束2,限制每个非零元素的取值范围为整数集P={0,±1,±2}。
通过枚举每一个对角矩阵中非零元素位置的所有组合,以及非零元素取值的所有组合,寻找使逼近误差最小的最优方案。
计算逼近误差时,采用矩阵范数比较DFT矩阵和分解矩阵乘积之间的差值。
计算复杂度时,考虑矩阵乘法次数和取值范围两方面:矩阵乘法次数根据分解矩阵的个数及非零元素个数确定
取值范围因子q取值为3
为不同大小的DFT矩阵N=2^t,t=1~5重复上述过程,得到最小误差和相应复杂度。
设DFT矩阵为$F_N$,要将其逼近为K个对角矩阵$D_k$的乘积:
$$F_N \approx D_1D_2...D_K$$
其中每个$D_k$满足:
- 非零元素数量 not more than 2 (约束条件1)
- 非零元素取值范围为整数集P(约束条件2)
则逼近过程为:
(1) 枚举$D_k$中非零元素位置的所有组合:
$pos_k = (i, j), i\neq j, i,j=1,...,N$
(2) 对每个组合,枚举非零元素的取值范围:
$D_k[i,i] \in P, D_k[j,j] \in P$
(3) 计算每个取值组合下的逼近误差:
$error = \frac{1}{N}|F_N - D_1D_2...D_K|_F$
(4) 选择使error最小的非零元素位置和取值组合
(5) 计算复杂度:
$C = q \times L$
其中$q=3$是取值范围因子,$L$为矩阵乘法次数
| import numpy as np from itertools import combinations def dft_matrix(N): i, j = np.meshgrid(np.arange(N), np.arange(N)) omega = np.exp(-2 * np.pi * 1j / N) W = np.power(omega, i * j) return W / np.sqrt(N) def diagonal_matrix(N, pos, values): D = np.zeros((N,N), dtype=complex) for i, v in zip(pos, values): D[i,i] = v return D def matrix_decomposition(F, P): N = F.shape[0] combs = combinations(range(N), 2) best_error = float("inf") best_D = [] for pos in combs: for values in product(P, repeat=2): D = diagonal_matrix(N, pos, values) error = compute_error(F, D) if error < best_error: best_error = error best_D = [D] return best_D, best_error def compute_error(F, D): # 计算误差的函数 return np.linalg.norm(F - D, 'fro') / np.sqrt(F.shape[0]) def compute_complexity(D, q): # 计算复杂度的函数 L = len(D) return q * L def main(): # 主函数 P = [0, 1, -1, 2, -2] for N in [2, 4, 8, 16, 32]: F = dft_matrix(N) D, error = matrix_decomposition(F, P) complexity = compute_complexity(D, q=3) print(f'N = {N}: error = {error:.4f}, complexity = {complexity}') if __name__ == '__main__': main() |
问题4
研究对Kronecker积矩阵的低复杂度逼近。当N1=4,N2=8时,具体思路如下:
- 根据定义,Kronecker积矩阵可以表示为:
F_N = F_4 ⊗ F_8
- 分别对F_4和F_8进行适当的低秩矩阵分解:
F_4 ≈ D_1D_2...D_m
F_8 ≈ E_1E_2...E_n
- 然后根据Kronecker积的性质,有:
F_N ≈ (D_1D_2...D_m) ⊗ (E_1E_2...E_n)
= (D_1⊗E_1)(D_2⊗E_2)...(D_m⊗E_n)
- 矩阵D和E的分解要满足稀疏性约束和取值范围约束。
- 通过搜索找到使逼近误差最小的D和E的分解。
- 计算复杂度时考虑D、E中矩阵的数目及稀疏性。
设$F_N$为$N=N_1N_2$阶的Kronecker积矩阵:
$F_N=F_{N_1}\otimes F_{N_2}$
其中$F_{N_1}$和$F_{N_2}$分别是$N_1$阶和$N_2$阶DFT矩阵。
对$F_{N_1}$和$F_{N_2}$分别进行低秩分解:
$F_{N_1}\approx D_1D_2\cdots D_M$
$F_{N_2}\approx E_1E_2\cdots E_L$
其中矩阵$D_i,E_j$满足约束条件:
- 每行最多2个非零元素(约束1)
- 非零元素取值范围为整数集P(约束2)
则根据Kronecker积的性质,有:
$F_N\approx(D_1D_2\cdots D_M)\otimes(E_1E_2\cdots E_L)$
$=(D_1\otimes E_1)(D_2\otimes E_2)\cdots(D_M\otimes E_L)$
搜索找到使逼近误差最小的$D_i,E_j$的最优分解,然后计算相应的复杂度。
Kronecker积矩阵保留了被合成矩阵的结构特征,这为低秩逼近提供了可能。
将大维DFT矩阵分解为多个小维DFT矩阵的Kronecker积,可以分别对小矩阵进行逼近,降低了优化难度。
将逼近问题分解为多个小规模子问题,符合“分治”的一般思想。
Kronecker积运算保留了矩阵乘法,可以继续使用低秩矩阵分解逼近的思路。
分解出的小矩阵满足稀疏性约束,可以有效减少乘法复杂度。
小矩阵取值范围限制也降低了每个乘法的计算复杂度。
可以通过搜索找到最小逼近误差的小矩阵分解,保证一定的逼近精度。
矩阵分解数量和取值范围可根据精度需求调整,实现可配置化。
| import numpy as np from scipy.linalg import kron def dft_matrix(n): i, j = np.meshgrid(np.arange(n), np.arange(n)) omega = np.exp(-2 * np.pi * 1j / n) W = np.power(omega, i*j) return W / np.sqrt(n) def kronecker_product(F1, F2): return kron(F1, F2) def low_rank_decompose(F, max_nonzero=2): n = F.shape[0] D = [] for i in range(n): d = np.diag([F[i,i]] + [0]*(n-1)) D.append(d) D_comb = list(combinations(D, max_nonzero)) # 选择误差最小的组合 F_approx = np.identity(n) for d in D_comb[best_index]: F_approx = F_approx @ d error = np.linalg.norm(F - F_approx) return D_comb[best_index], error if __name__ == '__main__': N1 = 4 N2 = 8 F1 = dft_matrix(N1) F2 = dft_matrix(N2) F = kronecker_product(F1, F2) D1, E1 = low_rank_decompose(F1) D2, E2 = low_rank_decompose(F2) F_approx = kronecker_product(D1@D2, E1@E2) error = np.linalg.norm(F - F_approx) / (N1*N2) print(error) |
问题五、
增加了精度限制要求,RMSE≤0.1。这个问题的主要难点是需要在满足精度约束的前提下,通过调整矩阵分解中元素的取值范围,来获得最小的硬件复杂度。针对这个问题,具体思路是:
- 使用问题3中矩阵分解的方法,将DFT矩阵F_N分解为多个对角矩阵的乘积。
- 对取值范围P进行递增搜索,比如依次取[0,±1]、[0,±1,±2]等,直到满足精度要求。
- 在每个取值范围下,搜索非零元素位置和取值,使RMSE最小。
- 记录下满足精度要求的最小取值范围。
- 在这个取值范围下,计算相应的硬件复杂度。
- 对不同大小的DFT矩阵N重复上述过程。
设DFT矩阵为$F_N$,将其分解为K个对角矩阵$D_k$的乘积:
$$F_N \approx D_1D_2\cdots D_K$$
其中每个$D_k$满足:
- 每行最多2个非零元素(约束1)
- 非零元素取值范围为整数集$P$(约束2)
要使逼近误差满足要求:
$$\text{RMSE} = \frac{|F_N - D_1D_2\cdots D_K|_F}{N} \leq 0.1$$
进行以下迭代搜索:
- 初始化取值范围:$P={0, ±1, ±2}$
- 在当前$P$下,搜索$D_k$的最优分解,使RMSE最小
- 如果RMSE $> 0.1$,扩大取值范围$P$,增加整数集大小
- 重复2)3),直到RMSE $\leq0.1$
- 输出此时的$P$和对应的复杂度$C$
其中,复杂度计算如前。
通过调整取值范围,可以满足精度要求,并使复杂度尽可能小。
设置精度约束RMSE≤0.1是问题的实际需求,方法必须首先满足这一约束。
通过搜索逐步扩大取值范围,可以系统地满足精度需求。
取值范围最小时,对应复杂度也最小,所以可以找到复杂度最小的解。
矩阵分解方法满足稀疏性,可以减少乘法次数,降低复杂度。
取值范围小,可以减少单个乘法的计算量,也降低了复杂度。
搜索可以找到精度和复杂度的最优trade-off。
不同大小矩阵可统一适用该方法,具有普适性。
可以获得在给定精度需求下的最小复杂度方案。
矩阵分解个数、取值范围都可配置,实现灵活可控。
| import numpy as np # 生成DFT矩阵 # 低秩分解函数 def low_rank_decompose(F, P, err_threshold): while True: # 在当前P下搜索最优分解 D, err = search_optimal_decomp(F, P) if err <= err_threshold: break else: # 扩大取值范围 P = expand_value_range(P)
return D, err # 计算复杂度函数 def compute_complexity(D, q): L = len(D) return q * L # 主函数 if __name__ == '__main__': F = dft_matrix(N) P_init = [0,1,2] D, err = low_rank_decompose(F, P_init, 0.1) q = len(P) comp = compute_complexity(D, q) |
消融实验分析:
基准模型:使用完整的方法,即矩阵分解+取值范围控制+稀疏性约束。测量其逼近误差RMSE和复杂度C。
移除取值范围控制:仅使用矩阵分解+稀疏性约束,不限制取值范围。测量RMSE和C。
移除稀疏性约束:仅使用矩阵分解+取值范围控制,不要求稀疏性。测量RMSE和C。
仅矩阵分解:不使用取值范围控制和稀疏性约束。测量RMSE和C。
对比不同模型的RMSE和C。高RMSE表示逼近精度损失;高C表示复杂度增加。
矩阵分解是实现低复杂度逼近DFT的有效方法,但需要设计实现稀疏性。
约束矩阵中元素的取值范围,可以降低单个乘法的计算量。
在满足精度需求前提下,通过搜索可以找到使复杂度最小的分解方案。
对Kronecker积矩阵进行分解,可以将大型DFT分解为多个小矩阵,降低优化难度。
消融实验可以验证不同设计决策对逼近误差和复杂度的影响。
需要权衡误差精度与计算复杂度,根据实际需求确定可接受的trade-off。
该方法可以作为一种替代FFT的低复杂度DFT实现策略。
优化搜索和代码实现等细节亟待进一步改进。
