华为杯数学建模B题

当大家面临着复杂的数学建模问题时，你是否曾经感到茫然无措？作为2021年美国大学生数学建模比赛的O奖得主，我为大家提供了一套优秀的解题思路，让你轻松应对各种难题。
让我们来看看研赛的B题呀~！

问题重述

DFT在通信等领域的重要应用,以及目前采用FFT计算DFT的硬件开销大的问题。提出了将DFT矩阵分解为整数矩阵乘积逼近的方法来降低硬件复杂度。
建模目标是对给定的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=15$给出不同约束条件下的优化问题,要求求出最小RMSE和相应的硬件复杂度C。

问题一：

要求在约束条件1(每个矩阵最多2个非零元素)下,对DFT矩阵$F_N(N=2^t,t=1,2,\mathrm{3...})$进行分解逼近,并计算最小误差和硬件复杂度。
这里采用的思路是:
1.将DFT矩阵F_N拆分为多个对角矩阵的乘积,每个对角矩阵只有一个非零元素,这样就满足了约束条件1。
2.对角矩阵的顺序和元素值可以通过搜索算法优化,以得到最小的逼近误差。
3.由于本题中没有限制取值范围,为简化计算,可将所有非零元素设为1。
4.硬件复杂度即为矩阵乘法次数,这里每个矩阵只有一个非零元素,所以复杂度就是矩阵个数。
例如当N=4时:
$F_4\approx$ $\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 0\end{array}\right]\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 1\end{array}\right]$

按此方法,计算了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线性相关。

1.DFT矩阵F_N的定义:
$$F_N=\frac{1}{\sqrt{N}}\left[\begin{array}{ccccccccccccccccc}1& 1& 1& \cdots & 11& w& w^2& \cdots & w^{N-1}⋮& ⋮& ⋮& \ddots & ⋮1& w^{N-1}& w^{2(N-1)}& \cdots & w^{(N-1)(N-1)}\end{array}\right]$$其中$w=e^{-j2\pi /N}$。
2.将F_N拆分为N个对角矩阵的乘积:
$$F_N\approx D_1{D}_2\cdots D_N$$
其中$D_k$为仅第k个对角元素为1的对角矩阵:
$$D_k=\left[\begin{array}{ccccccccccccccccc}0& & & \ddots & & & 1_{kk}& & & & & \ddots & & & & & 0\end{array}\right]$$
3.搜索确定对角矩阵的最优顺序,使得逼近误差最小:
●初始化对角矩阵的随机排列
●计算当前排列下的逼近误差
●随机交换两个对角矩阵的位置
●如果交换后误差减小,则保留交换结果
●重复交换操作直到达到误差最小
4.逼近误差的计算:
$$RMSE=\frac{1}{N}\sqrt{|F_N-D_1{D}_2\cdots D_N{|}_F^2}$$
5.硬件复杂度即为矩阵乘法次数,这里每个D_k矩阵仅有一个非零元素,所以复杂度就是矩阵个数N。
6.按此方法,计算从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。

1.F_4 的定义如下:
$$F_4=\frac{1}{2}\left[\begin{array}{ccccccccccccc}1& 1& 1& 11& j& -1& -j1& -1& 1& -11& -j& -1& j\end{array}\right]$$
2.将其分解为4个对角矩阵Di:
$$F_4\approx D_1{D}_2{D}_3{D}_4$$
其中Di是仅第i个对角元素非零的对角矩阵。
3.根据元素取值范围P={0,±1,±2},对Di的非零元素取值进行穷举,选择误差最小的取值:
$$\begin{array}{rl}{D}_1& =\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\\ D_2& =\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\\ D_3& =\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 0\end{array}\right]\\ D_4& =\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 1\end{array}\right]\end{array}$$
4.逼近误差计算:
$$RMSE=\frac{1}{4}|\frac{1}{2}{F}_4-D_1{D}_2{D}_3{D}_4{|}_F=\frac{1}{2}$$
5.计算复杂度:
●每个矩阵乘法都包含一个复数乘法。
●根据元素取值范围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_1{D}_2...D_K$$
其中每个$D_k$满足:
1.非零元素数量 not more than 2 (约束条件1)
2.非零元素取值范围为整数集P(约束条件2)
则逼近过程为:
(1) 枚举$D_k$中非零元素位置的所有组合:
$po{s}_k=(i,j),i\ne 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_1{D}_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时,具体思路如下:
1.根据定义,Kronecker积矩阵可以表示为:
$$F_N=F_4\otimes F_8$$
2.分别对F_4和F_8进行适当的低秩矩阵分解:
$$\begin{gathered} F_4\approx D_1{D}_2...D_m \\ {F}_8\approx E_1{E}_2...E_n \end{gathered}$$
3.然后根据Kronecker积的性质,有:
$$\begin{gathered} F_N\approx (D_1{D}_2...D_m)\otimes (E_1{E}_2...E_n) \\ =(D_1\otimes E_1)(D_2\otimes E_2)...(D_m\otimes E_n) \end{gathered}$$
4.矩阵D和E的分解要满足稀疏性约束和取值范围约束。
5.通过搜索找到使逼近误差最小的D和E的分解。
6.计算复杂度时考虑D、E中矩阵的数目及稀疏性。

设$F_N$为$N=N_1{N}_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_1{D}_2\cdots D_M$
$F_{N_2}\approx E_1{E}_2\cdots E_L$
其中矩阵$D_i,E_j$满足约束条件:
1.每行最多2个非零元素(约束1)
2.非零元素取值范围为整数集P(约束2)
则根据Kronecker积的性质,有:
$F_N\approx (D_1{D}_2\cdots D_M)\otimes (E_1{E}_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。这个问题的主要难点是需要在满足精度约束的前提下,通过调整矩阵分解中元素的取值范围,来获得最小的硬件复杂度。针对这个问题,具体思路是:
1.使用问题3中矩阵分解的方法,将DFT矩阵F_N分解为多个对角矩阵的乘积。
2.对取值范围P进行递增搜索,比如依次取[0,±1]、[0,±1,±2]等,直到满足精度要求。
3.在每个取值范围下,搜索非零元素位置和取值,使RMSE最小。
4.记录下满足精度要求的最小取值范围。
5.在这个取值范围下,计算相应的硬件复杂度。
6.对不同大小的DFT矩阵N重复上述过程。

设DFT矩阵为$F_N$,将其分解为K个对角矩阵$D_k$的乘积:
$$F_N\approx D_1{D}_2\cdots D_K$$
其中每个$D_k$满足:
1.每行最多2个非零元素(约束1)
2.非零元素取值范围为整数集$P$(约束2)
要使逼近误差满足要求:
$$\text{RMSE}=\frac{|F_N-D_1{D}_2\cdots D_K{|}_F}{N}\le 0.1$$
进行以下迭代搜索:
1.初始化取值范围:$P=0,\pm 1,\pm 2$
2.在当前$P$下,搜索$D_k$的最优分解,使RMSE最小
3.如果RMSE $>0.1$,扩大取值范围$P$,增加整数集大小
4.重复2)3),直到RMSE $\le 0.1$
5.输出此时的$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实现策略。
优化搜索和代码实现等细节亟待进一步改进。

想要获取更多完整版看看这里哈~
(5 封私信 / 2 条消息) 如何评价2023数学建模研赛B题？ - csdn