2023年全国研究生数学建模竞赛华为杯

B题 DFT类矩阵的整数分解逼近

原题再现：

  一、问题背景
  离散傅里叶变换（Discrete Fourier Transform，DFT）作为一种基本工具广泛应用于工程、科学以及数学领域。例如，通信信号处理中，常用DFT实现信号的正交频分复用（Orthogonal Frequency Division Multiplexing，OFDM）系统的时频域变换（见图1）。另外在信道估计中，也需要用到逆DFT（IDFT）和DFT以便对信道估计结果进行时域降噪（见图2）。

  在芯片设计中，DFT计算的硬件复杂度与其算法复杂度和数据元素取值范围相关。算法复杂度越高、数据取值范围越大，其硬件复杂度就越大。目前在实际产品中，一般采用快速傅里叶变换（Fast Fourier Transform，FFT）算法来快速实现DFT，其利用DFT变换的各种性质，可以大幅降低DFT的计算复杂度（参见[1][2]）。然而，随着无线通信技术的演进，天线阵面越来越大，通道数越来越多，通信带宽越来越大，对FFT的需求也越来越大，从而导致专用芯片上实现FFT的硬件开销也越大。为进一步降低芯片资源开销，一种可行的思路是将DFT矩阵分解成整数矩阵连乘的形式。



  可以看到在该方案中，分解后的矩阵元素均为整数，从而降低了每个乘法器的复杂度；另外A_1~A_4的稀疏特性可以减少乘法运算数量。可以看出，这其实是一种精度与硬件复杂度的折中方案，即损失了一定的计算精度，但是大幅降低了硬件复杂度。在对输出信噪比要求不高的情况下可以优先考虑此类方案。

  目前使用FFT进行DFT计算的方案硬件复杂度较高，因为我们希望研究一种替代方案来降低DFT计算的硬件复杂度，但同时我们对精度也有一定要求。请针对以下问题分别设计分解方法，既能最小化RMSE，同时又使得乘法器的数量尽量少。
  A中矩阵的个数K的取值并没有限制，也是优化的变量之一。但需要注意，一般情况下，K越小，硬件复杂度越低，但是如果增加矩阵的个数可以使得矩阵中包含更多的简单元素（0、±1、±j或(±1±j)），硬件复杂度也可能会降低，因此，需要根据硬件复杂度C的定义合理的设计K。

  问题1：首先通过减少乘法器个数来降低硬件复杂度。由于仅在非零元素相乘时需要使用乘法器，若A_k矩阵中大部分元素均为0，则可减少乘法器的个数，因此希望A_k为稀疏矩阵。对于N=2^t,t=1,2,3,…的DFT矩阵F_N，请在满足约束1的条件下，对最优化问题(6)中的变量A和β进行优化，并计算最小误差（ 即(6)的目标函数，下同）和方案的硬件复杂度C（由于本题中没有限定A_k元素的取值范围，因此在计算硬件复杂度时可默认q=16）。

  问题2：讨论通过限制A_k中元素实部和虚部取值范围的方式来减少硬件复杂度的方案。对于N=2^t,t=1,2,3,4,5的DFT矩阵F_N，请在满足约束2的条件下，对A和β进行优化，并计算最小误差和方案的硬件复杂度C。

  问题3：同时限制A_k的稀疏性和取值范围。对于N=2^t,t=1,2,3,4,5的DFT矩阵F_N，请在同时满足约束1和2的条件下，对A和β进行优化，并计算最小误差和方案的硬件复杂度C。

  问题4：进一步研究对其它矩阵的分解方案。考虑矩阵F_N=F_N1⊗F_N2，其中F_N1 和F_N2分别是N_1和N_2维的DFT矩阵，⊗表示Kronecker积（注意F_N非DFT矩阵）。当N_1=4, N_2=8时，请在同时满足约束1和2的条件下，对A和β进行优化，并计算最小误差和方案的硬件复杂度C。

  问题5：在问题3的基础上加上精度的限制来研究矩阵分解方案。要求将精度限制在0.1以内，即RMSE≤0.1。对于N=2^t,t=1,2,3…的DFT矩阵F_N，请在同时满足约束1和2的条件下，对A和β,P进行优化，并计算方案的硬件复杂度C。

  附录一：名词解释
  复数乘法次数/复乘次数：进行复数乘法的次数，例如(1+2j)×(2+2j)为一次复乘。
  硬件复杂度：本题中，仅考虑乘法器带来的硬件复杂度，硬件复杂度仅与乘法器个数和每个乘法器的复杂度相关
  乘法器个数：本题中，乘法器个数即为复乘次数
  单个乘法器的复杂度：单个乘法器的复杂度与乘法器的设计方法和输入数据的位宽等因素相关。在本题中，将乘法器的复杂度简化为仅与输入数据的取值范围相关。对于复数g∈{x+jy│x,y∈P},P={0,±1,±2,…,±2^(q-1) },其与任意复数z相乘的复杂度为q。

整体求解过程概述(摘要)

  离散傅里叶变换（Discrete Fourier Transform，DFT）傅里叶分析方法是信号分析的最基本方法，傅里叶变换是傅里叶分析的核心，通过它把信号从时间域变换到频率域，进而研究信号的频谱结构和变化规律。在芯片设计中，DFT计算的硬件复杂度与其算法复杂度和数据元素取值范围相关。算法复杂度越高、数据取值范围越大，其硬件复杂度就越大。常规的降低硬件复杂度的方法如快速傅里叶变换已经渐渐无法满足芯片日益增长的需求。因此，急需设计新的算法，此算法不仅能是误差控制在一定范围内，又能有效降低 DFT 过程带来的硬件复杂度过大的问题。
  针对第一问，本问不用考虑分解后矩阵A的取值范围的问题，只需要满足分解后的矩阵每行至多只有2个非零元素，以及最小误差尽可能低的约束条件下，计算出近似矩阵的硬件复杂度。本问选择了三种模型来计算 DFT 矩阵的最小误差和硬件复杂度。奇异值分解法通过特征向量将DFT矩阵进行分解，同时还达到了降维的目的，通过将DFT矩阵分解成三个矩阵来使误差达到最小；分块矩阵分解法通过将 DFT 矩阵分解成一个个小的分块矩阵，因为维数越小分块子矩阵形式越简单，误差显著降低；矩阵乘法拟合则利用 DFT 矩阵的对称性和穷举法将矩阵分解成多个矩阵连乘的形式。
  针对第二问，本问不用考虑每行非零元素个数的问题，只需要满足分解后的矩阵每个元素的取值范围，以及在N=2，4，8，16，32的情况以及最小误差尽可能低的约束条件下，计算出近似矩阵的硬件复杂度。本问选择了三种模型来计算 DFT 矩阵的最小误差和硬件复杂度。蝶形运算分解法运用了快速傅里叶变换的思想，在此基础上利用蝶形变换将DFT矩阵进行分解；分块矩阵分解法和矩阵乘法拟合在解决问题二是仍然适用。
  针对第三问，本问需要同时考虑每行非零元素个数的以及分解后的矩阵每个元素的取值范围的约束条件，在此基础上，在最小误差尽可能低的约束条件下，计算出近似矩阵的硬件复杂度。本问选择了两种模型来计算 DFT 矩阵的最小误差和硬件复杂度，即分块矩阵分解法和矩阵乘法拟合。这两种方法在多种约束条件下仍能够很好的解决问题。
  针对第四问，本问需要考虑如何对4点DFT矩阵与8点DFT矩阵的Kronecker积的矩阵FN进行矩阵分解，在最小误差尽可能低的约束条件下，计算出近似矩阵的硬件复杂度。本问选择通过SVD+穷举法的模型来计算DFT矩阵的最小误差和硬件复杂度。首先对4点DFT矩阵与8点DFT矩阵进行SVD分解，之后利用Kronecker积的性质将FN矩阵转化为多个矩阵相乘的形式。
  针对第五问，本问需要同时考虑每行非零元素个数的以及分解后的矩阵每个元素的取值范围的约束条件，在此基础上增加将精度限制在0.1以内，即RMSE≤0.1的要求，计算出近似矩阵的硬件复杂度。本问选择分块矩阵分解模型对第三问结果矩阵进行进一步优化，在满足题目要求下，计算出分解后矩阵的最小误差和硬件复杂度。

模型假设：

  1. 硬件复杂度仅考虑乘法器的复杂度，硬件复杂度与乘法器个数和每个乘法器的复杂度成正比。
  2. 单个乘法器的复杂度简化为仅与复数中实部和虚部的取值范围相关。
  3. 乘法器个数定义为复数乘法的次数，且与0,±1,±j±1±j相乘时不计入乘法次数。

问题分析：

  DFT 是对信号向量进行线性正交变换的一个过程，对一个 N 维的信号直接进行 N 点 DFT 需要进行N2次复数乘法和 N(N − 1)次加法，对于硬件复杂度和算力要求较高。所以需 要设计快速算法将 DFT 计算成本降低，通常使用的方式为使用快速傅里叶变换算法（FFT） 来实现 DFT 过程但 FFT 的硬件开销不满足现有需求需要设计方案降低 FFT 的硬件开销。 一种合理的方式为使用矩阵连乘来拟合 DFT 矩阵，通过设计合理的矩阵可以降低 FFT 中 的硬件开销并满足精度要求。 FFT 中的硬件复杂度主要由矩阵连乘中复数乘法的次数和乘法器的复杂度决定，减少 复数乘法的次数并降低每个乘法器的复杂度可以有效的使得 FFT 中的硬件开销变小。连乘 中矩阵尽可能稀疏会使得乘法的次数有效减少，而乘法器的复杂度可以通过减少乘法元素 的有效位数有效的降低乘法器的复杂度。因此为降低 FFT 中的硬件开销，需要对连乘中矩 阵的稀疏性和元素有效位数进行合理设计，在精度尽可能大的条件下硬件复杂度降低。
  在设计好的矩阵稀疏性和元素有效位数约束下，根据不同维度的 DFT 矩阵来优化需要 连乘的矩阵个数、每个矩阵元素的值和矩阵缩放因子来得到最优的矩阵连乘形式。但在规定矩阵元素虚部和实部均为整数的条件下，直接对每个矩阵元素中的值进行优化为一个庞 大的多元整数优化问题，不是一个明智的思路。因此需要分析 DFT 矩阵的性质和旋转因子Wn的性质并利用这些性质设计合理的优化方案，从而在可能的范围内求解出最合适的矩阵连乘形式。

模型的建立与求解整体论文缩略图

全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可

程序代码：

部分程序代码如下:

NNN=8;
 O_FN=zeros(1,NNN);
 O_FN_decom= zeros(1,NNN);
 for i=3:NNN+2
 N=2^i;
 FN=produce_DFT(N);
 [FN_decom,N_length]=my_CT(N);
 O_FN(i−2)= calculate_On(FN,N);
 for j=1:N_length
 O_FN_decom(i−2)=O_FN_decom(i−2)+calculate_On(FN_decom{j},
 N);
 end
 end
 semilogy(2.^(3:NNN+2),O_FN,'−*')
 gridon
 holdon
 semilogy(2.^(3:NNN+2),O_FN_decom,'−*')
 for i=1:NNN
 line([2^(i+2),2^(i+2)],[O_FN(i),O_FN_decom(i)],'LineStyle','−−'
 ,'Color','g');
 end
 xlabel('N:矩阵维数')
 ylabel('C:硬件复杂度')
 legend('C(F_N)','C(F_N的精确分解)')

function result= produce_DFT(N)
 result=zeros(N,N);
 w=exp(((−2*pi)/N)*1i);
 result(:,1)=ones(N,1);
 temp=ones(N,1);
 for i=2:N
 temp(i)=w*temp(i−1);
 end
 result(:,2)=temp;
 for i=3:N
 result(:,i)=result(:,i−1).*temp;
 end

function [result,t]=my_CT(N)
 w=exp(((−2*pi)/N)*1i);
 t=log2(N);
 result=cell(1,t+1);
 P_eye= eye(N);
 for i=1:t
 temp= zeros(N,N);
 P_temp= zeros(N,N);
 Ni=N/(2^(i−1));
 Ni_2= Ni/2;
 re_temp=zeros(Ni,Ni);
 P=zeros(Ni,Ni);
 for j=1:Ni/2
 re_temp(j,j)=1;
 re_temp(j,j+Ni_2)= w^(2^(i−1)*(j−1));
 P(j,2*j−1)=1;
 end
 for j=Ni_2+1:Ni
 re_temp(j,j−Ni_2)= 1;
 re_temp(j,j)=−w^(2^(i−1)*(j−1−Ni_2));
 P(j,2*(j−Ni_2))=1;
 end
 for j=1:2^(i−1)
 temp(1+Ni*(j−1):Ni*j,1+Ni*(j−1):Ni*j)=re_temp;
 P_temp(1+Ni*(j−1):Ni*j,1+Ni*(j−1):Ni*j)=P;
 end
 result{i}= temp;
 P_eye= P_temp*P_eye;
 end
 result{t+1}=P_eye;
 t=t+1;

全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可