2023 “华为杯”研赛B题，建模解析，小鹿学长带队指引全代码文章与思路

小鹿学长联合20多位国二以上的好队友，共同提供了进阶版精华包呀，有兴趣的同学可以点击下面链接看看～～
link：（进阶版）2023 研赛B题，华为杯建模解析，小鹿学长带队指引全代码文章与思路

问题重述
本题研究如何用多个整数矩阵的乘积来逼近DFT矩阵,目的是用这种方法代替目前在芯片上用于DFT计算的FFT算法,以降低硬件复杂度。

给定N维DFT矩阵F_N,要求设计K个矩阵A_1到A_K,使得它们的乘积最接近于βF_N,其中β是一个缩放系数。目标是最小化它们之间的Frobenius范数误差。

硬件复杂度C定义为:乘法器个数q × 复数乘法次数L。其中q表示矩阵元素的取值范围,L表示进行复数乘法的次数。

问题要求:

在矩阵行数限制为2的条件下,优化A和β以最小化误差,计算最小误差和硬件复杂度C。
在矩阵元素取值范围限制为整数的条件下,优化A和β,计算误差和C。
在同时考虑稀疏性约束和取值范围约束的条件下,优化A和β,计算误差和C。
对Kronecker矩阵乘积情况进行分析。
在一定精度限制下,优化A、β和取值范围q,计算C。

目标是在保证一定的计算精度的前提下,尽可能减少硬件复杂度C。

问题一
问题1要求在约束每个A_(k)矩阵的每行最多有2个非零元素的条件下,优化变量A和β以最小化DFT矩阵F_N与A_1A_2…A_K乘积之间的Frobenius范数误差。

建模思路:

定义优化变量:

A = {A_1,A_2,…,A_K},其中A_k是N阶方阵
β,实数缩放因子

定义优化目标:

\min_{\mathcal{A},\beta}\frac{1}{N}\sqrt{\|\beta F_N - A_1A_2\cdots A_K\|_F^2}
minA,β

定义约束条件:

对于k=1,2,…,K:

A_k的每行有至多2个非零元素

计算复杂度:

默认q=16

则C = q\times L

C=q×L

其中L为进行复数乘法的次数,与0、±1、±j或(±1±j)相乘不计入L。

求解优化问题,得到最优A、β,计算最小误差及复杂度C。

import numpy as np
from scipy.optimize import minimize

# DFT矩阵定义
N = 8
F_N = np.zeros((N, N), dtype=complex)
for k in range(N):
for n in range(N):
F_N[k, n] = np.exp(-2j * np.pi * k * n / N) / np.sqrt(N)

# 优化目标函数
def optimize_func(x, F_N):
A = x[:-1].reshape(K, N, N) # 构造A矩阵
beta = x[-1]
F_approx = beta * F_N
for i in range(K):
F_approx = np.dot(F_approx, A[i]) # 逐个矩阵乘积
loss = np.linalg.norm(F_N - F_approx, 'fro') / N
return loss

# 约束条件函数
def constraint_func(x):
A = x[:-1].reshape(K, N, N)
for i in range(K):
if np.count_nonzero(A[i], axis=1).max() > 2:
return False
return True

# 初始化参数
K = 3
init_x = [np.random.randn(N, N) for i in range(K)]
问题二

定义优化变量:

A = \{A_1,A_2,\ldots,A_K\}

A={A1
​,A2
​,…,AK
​},其中A_k\in \mathbb{C}^{N\times N}

Ak
​∈CN×N
\beta\in \mathbb{R}

β∈R,缩放因子

定义优化目标:

\min_{\mathcal{A},\beta}\frac{1}{N}\|\beta F_N - A_1A_2\cdots A_K\|_F

定义约束条件:

对于k=1,2,\ldots,K

k=1,2,…,K:

\begin{equation}
A_k[i,j] \in {x+jy|x,y\in P}, \quad i,j=1,2,\ldots,N
\end{equation}
其中P={0,±1,±2,…,±2q−1},此处取q=3
q=3

计算复杂度

(1) 定义L_k为矩阵A_k进行复数乘法的次数。

(2) L_k的计算:遍历矩阵A_k​的每个元素A_k[i,j],判断该元素是否属于{0, ±1, ±j, (±1±j)}这些简单乘法的元素,如果不是,则L_k加1。

(3) 重复上述操作,可以得到每个矩阵A_k​的复数乘法次数L_kLk​。

(4) 根据复杂度公式C = q\times(\sum_{k=1}^{K} L_k)
C=q×(∑k=1
K
​Lk
​),可以求出总的复杂度C

求解优化问题

(1) 将优化目标函数和约束条件建立成一个无约束优化问题。
(2) 初始化优化变量A和β。
(3) 使用优化算法(如梯度下降法、牛顿法等)求解这个无约束优化问题。
(4) 比较不同初始化的优化结果,选择其中误差最小的解作为最优解。
(5) 将最优解A,\beta
A,β代入公式,计算最小误差和复杂度C

import numpy as np
from scipy.optimize import minimize

N = 8
F_N = np.fft.fft(np.eye(N)) / np.sqrt(N)

# 优化目标函数
def optimize_func(x, F_N):
A = x[:-1].reshape(K, N, N)
beta = x[-1]
loss = np.linalg.norm(beta*F_N - np.prod(A, axis=0), 'fro') / N
return loss

# 约束条件函数
def constraint_func(x):
A = x[:-1].reshape(K, N, N)
for i in range(K):
if np.any(np.abs(A[i]) > 2**2):
return False
return True