\section{Introduction}

Federated learning (FL) is a promising technique to reduce the communication overhead while protecting the data privacy at users for effective machine learning \cite{fed1}.
Although FL was originally developed for wire-line connected systems \cite{fed1}, the development of mobile intelligent systems such as autonomous driving \cite{zijian} calls for wireless connections between the server and users, giving rise to a new paradigm termed wireless FL or edge FL \cite{fed2,fed3,fed4,air2,air3,air4,air5}.
However, the convergence of wireless FL may require an exceedingly long time due to limited capacity of wireless channels during the uplink model aggregation step.
To reduce the transmission delay, various wireless FL designs have been proposed, which are mainly categorized into digital modulation \cite{fed2,fed3,fed4} and analog modulation \cite{air2,air3,air4,air5} methods.

For digital modulation systems, data from different users are multiplexed either in the time or the frequency domain \cite{massive}.
Current works on delay reduction focus on reducing 1) the number of model aggregation iterations \cite{fed2}, 2) the number of users \cite{fed3}, or 3) the number of bits for representing the gradient of back propagation in each iteration \cite{fed4}.
However, since these strategies involve various approximations of the FL procedure, the performance of learning would be degraded inevitably.
On the other hand, the key advantage of analog modulation \cite{air2,air3,air4,air5} over digital modulation arises from the ground-breaking idea of over-the-air computation (AirComp).
Specifically, if multiple users upload their local parameters simultaneously, a superimposed signal, which represents a weighted sum of individual model parameters, is observed at the edge server.
By performing the minimum mean square error (MMSE) detection on the superimposed signal, an estimate of the global parameter vector can be obtained, thereby significantly shortening the required transmission time.
However, due to channel fading and noise in wireless systems, AirComp employed in single-antenna systems \cite{air2,air3} could result in large error in the estimation of global model parameters at the edge server, leading to slow convergence of FL iterations.
As a remedy, adopting MIMO beamforming \cite{air4,air5} could reduce the parameter transmission error by aligning the beams carrying the local parameters' information to the same desired spatial direction.
Nonetheless, the current transmit and receive beamforming designs in MIMO AirComp systems involve exceedingly high radio frequency (RF) chain costs and high computational complexities \cite{air4,air5}, preventing their practical implementation.

To fill the research gap, this paper proposes the unit-modulus wireless FL (UMWFL) framework which integrates AirComp and unit-modulus analog beamforming.
Specifically, the edge users possess a number of training data for local model updates.
Then, the trained model parameters are uploaded to the server via analog modulation.
To reduce the implementation costs, the edge server does not process the received model parameters at the baseband.
Instead, it applies a phase shift network in the RF domain and an RF sampling module to connect received antennas and transmit antennas for global model updates.
Upon receiving the broadcast, all users feed the received signals to an analog demodulator for parameter extraction.
Note that our UMWFL framework can significantly reduce the required RF chains in MIMO FL systems, thereby reducing the associated hardware and energy costs.
The proposed UMWFL framework also alleviates the straggler effects in both model uploading and broadcasting procedures.
On the other hand, despite the UMWFL problem being highly nonconvex, a large-scale optimization algorithm, termed penalty alternating minimization (PAM), is developed.
Experimental results show that the learning performance of the proposed PAM significantly outperforms the benchmark scheme without phase optimization.

\emph{Notation}: Italic letters, lowercase and uppercase bold letters represent scalars, vectors, and matrices, respectively.
Curlicue letters stand for sets and $|\cdot|$ is the cardinality of a set.
The operators $\|\cdot\|_2,(\cdot)^{T},(\cdot)^{H},(\cdot)^{-1}$ are the $\ell_2$-norm, transpose, Hermitian, and inverse of a matrix, respectively.
The operators $\partial f$ and $\nabla f$ are the partial derivative and the gradient of the function $f$.
The operators $\mathrm{vec}$ and $\mathrm{mat}$ denote the vectorization of a matrix and the matricization of a vector.
The function $[x]^+=\mathrm{max}(x,0)$, $\mathrm{Re}(x)$ takes the real part of $x$, $\mathrm{Im}(x)$ takes the imaginary part of $x$, $\mathrm{conj}(x)$ takes the conjugate of $x$, and $|x|$ is the modulus of $x$.
$\mathbf{I}_{N}$ denotes the $N\times N$ identity matrix and $\otimes$ denotes the Kronecker product.
Finally, $\mathrm{j}=\sqrt{-1}$, $\mathbb{E}(\cdot)$ denotes the expectation of a random variable and $\mathcal{O}(\cdot)$ is the big-O notation standing for the order of arithmetic operations.

\begin{figure}[!t]
 \centering
\includegraphics[width=0.85\textwidth]{system1.eps}
  \caption{The UMWFL system with a phase shifting edge server.}
\end{figure}

\section{Unit Modulus Wireless Federated Learning}

\setcounter{secnumdepth}{4}
We consider a UMWFL system shown in Fig.~1, which consists of an edge server equipped with $N$ antennas and $K$ single-antenna mobile users. The dataset and model parameter vector at user $k$ are denoted as $\mathcal{D}_k$ and $\mathbf{x}_k\in\mathbb{R}^{M\times 1}$, respectively.
Mathematically, the UMWFL procedure aims to solve the following optimization problem:
\begin{align}\label{FL}
\mathop{\mathrm{min}}_{\substack{\{\mathbf{x}_k\},\bm{\theta}}}
\quad&\underbrace{\frac{1}{\sum_{k=1}^{K}|\mathcal{D}_k|}\sum_{k=1}^K\sum_{\mathbf{d}_{k,l}\in\mathcal{D}_k}\Theta(\mathbf{d}_{k,l}, \bm{\theta})}_{:=\Lambda(\bm{\theta})}
\nonumber\\
\mathrm{s.t.}\quad&\mathbf{x}_1=\cdots=\mathbf{x}_K= \bm{\theta},
\end{align}
where $\Theta(\mathbf{d}_{k,l}, \bm{\theta})$ is the loss function corresponding to a single sample $\mathbf{d}_{k,l}$ ($1\leq l \leq |\mathcal{D}_k|$) in $\mathcal{D}_k$ given parameter vector $\bm{\theta}$, while $\Lambda(\bm{\theta})$ denotes the global loss function to be minimized.
The training of UMWFL model parameters (i.e., solving \eqref{FL}) in the considered edge system is naturally a distributed and iterative procedure, where each iteration involves four steps:
1) updating the local parameter vectors $(\mathbf{x}_1,\cdots,\mathbf{x}_K)$ using $\{\mathcal{D}_1,\cdots,\mathcal{D}_K\}$ at users $(1,\cdots,K)$, respectively;
2) transforming the local parameters $(\mathbf{x}_1,\cdots,\mathbf{x}_K)$ into transmit symbols $(\mathbf{s}_1,\cdots,\mathbf{s}_K)$ via power and phase designs;
3) aggregating $(\mathbf{x}_1,\cdots,\mathbf{x}_k)$ in an analog manner at the edge server;
and 4) broadcasting the results to users.
The details for the $i$-th iteration of UMWFL are elaborated below.

\subsubsection{Local Updating}
In the first step, let $\mathbf{x}^{[i]}_k(0)\in\mathbb{R}^{M\times 1}$ be the local parameter vector at user $k$ at the beginning of the $i$-th iteration.
To update $\mathbf{x}^{[i]}_k(0)$, user $k$ minimizes the loss function $\frac{1}{|\mathcal{D}_k|}\sum_{\mathbf{d}_{k,l}\in\mathcal{D}_k}\Theta(\mathbf{d}_{k,l}, \mathbf{x}_k)$ via the gradient descent method as
\begin{align}\label{sgd}
&\mathbf{x}^{[i]}_k(\tau+1)=\mathbf{x}^{[i]}_k(\tau)-\frac{\varepsilon}{|\mathcal{D}_k|}\sum_{\mathbf{d}_{k,l}\in\mathcal{D}_k}\nabla_{\mathbf{x}}\Theta\left(\mathbf{d}_{k,l},\mathbf{x}^{[i]}_k(\tau)\right),
\end{align}
where $\varepsilon>0$ is the step-size and $\tau$ is from $0$ to $E-1$ with $E$ being the number of local updates.
Then, $\{\mathbf{x}^{[i]}_k(E)|\forall k\}$ from all users are uploaded to the edge server.

\subsubsection{Parameter Uploading}
In the second step, user $k$ encodes its local parameter vector $\mathbf{x}^{[i]}_k(E)$ into a complex vector $\mathbf{s}^{[i]}_k$.
Since the model parameters in deep learning are real-valued numbers, in order to reduce the transmission time, every two model parameters are transmitted as a complex number.  That is,
\begin{align}
\mathbf{s}_k^{[i]}=&
\frac{\overbrace{\sqrt{p^{[i]}_k}\,\mathrm{e}^{\rm{j}\phi_k^{[i]}}}^{:=t_k^{[i]}}}{\sqrt{2\eta^{[i]}}}\,\Big[x_{k,1}^{[i]}(E)+\mathrm{j}\,x_{k,2}^{[i]}(E),\cdots,x_{k,M-1}^{[i]}(E)+\mathrm{j}\,x_{k,M}^{[i]}(E)\Big]^T, \label{sk}
\end{align}
where $p_k^{[i]}\in\mathbb{R}$ and $\phi_k^{[i]}$ are the transmit power and phase at user $k$, $\mathbf{s}_k^{[i]}\in\mathbb{C}^{S\times 1}$ with $S=M/2$, $x_{k,m}^{[i]}(E)$ is the $m$-th element of $\mathbf{x}_k^{[i]}(E)$.
The scaling factor $\eta^{[i]}$ is $\eta^{[i]}=\frac{1}{K}\sum_{k=1}^K\overline{\eta}_k^{[i]}$ with $\overline{\eta}_k^{[i]}=\frac{1}{M}\|\mathbf{x}_k^{[i]}(E)\|_2^2$ such that the average power of $\mathbf{s}_k^{[i]}$ is $\frac{1}{S}\mathbb{E}[\|\mathbf{s}_k^{[i]}\|_2^2]=p_k^{[i]}$.
To facilitate the subsequent derivations, we define the transmit coefficient $\{t_k^{[i]}\}$ in \eqref{sk} and $\{p_k^{[i]},\varphi_k^{[i]}\}$ can be recovered from $\{t_k^{[i]}\}$.

\subsubsection{Parameter Aggregation}
In the third step, the received signal $\mathbf{R}^{[i]}\in\mathbb{C}^{N\times S}$ at the server is
\begin{align}\label{Ri}
&\mathbf{R}^{[i]}=\sum_{k=1}^K\mathbf{h}_k^{[i]}(\mathbf{s}_{k}^{[i]})^T+\mathbf{Z}^{[i]},
\end{align}
where $\mathbf{h}_k^{[i]}\in\mathbb{C}^{N\times 1}$ is the uplink channel vector from user $k$ to the server and $\mathbf{Z}^{[i]}\in\mathbb{C}^{N\times S}$ is the matrix of the additive white Gaussian noise with covariance matrix
$\mathbb{E}\left[\mathrm{vec}(\mathbf{Z}^{[i]})\mathrm{vec}(\mathbf{Z}^{[i]})^H\right]=\sigma^2_b\mathbf{I}_{NS}$, where $\sigma^2_b$ is the noise power at the server.
Upon receiving the superimposed signal, the server processes $\mathbf{R}^{[i]}$ using a phase shift network $\mathbf{F}^{[i]}\in\mathbb{C}^{N\times N}$.
The phase shifted signals are then fed to the RF sampling module \cite{rfsampling} for data caching.
As shown in Fig.~1, the matrix $\mathbf{F}^{[i]}$ requires unit-modulus constraints on all of its elements \cite{analog,analog2}, i.e., $\mathbf{F}^{[i]}\in\mathcal{F}$ where
\begin{align}
&\mathcal{F}=\{\mathbf{F}:|F_{l,l'}^{[i]}|=1,\quad\forall l,l'\}. \label{setF}
\end{align}

\subsubsection{Parameter Decoding}
Finally, the server broadcasts $\sqrt{\gamma}\,\mathbf{F}^{[i]}\mathbf{R}^{[i]}\in\mathbb{C}^{N\times S}$ to all the users, where $\gamma>0$ is the power scaling factor adopted at the edge server.
The received signal at user $k$ is
\begin{align}
&\left(\mathbf{y}^{[i]}_k\right)^T=\left(\mathbf{g}_k^{[i]}\right)^H
\sqrt{\gamma}\,\mathbf{F}^{[i]}\mathbf{R}^{[i]}+\left(\mathbf{n}^{[i]}_k\right)^T, \label{yk}
\end{align}
where $\mathbf{g}_k^{[i]}\in\mathbb{C}^{N\times 1}$ is the downlink channel vector from the server to user $k$ and $\mathbf{n}_k^{[i]}\in\mathbb{C}^{S\times 1}$ is the vector of the additive white Gaussian noise with covariance matrix $\sigma^2_k\mathbf{I}_S$, where $\sigma^2_k$ is the noise power at user $k$.
User $k$ applies a normalization coefficient $r^{[i]}_k\in\mathbb{C}$ to $\mathbf{y}^{[i]}_k$ as
\begin{align}
\mathbf{x}^{[i+1]}_k(0)
=&
\sqrt{2\eta^{[i]}}\,\Big[\mathrm{Re}(r_k^{[i]}y^{[i]}_{k,1}),\mathrm{Im}(r^{[i]}_ky^{[i]}_{k,1}),\cdots,\mathrm{Re}(r^{[i]}_ky^{[i]}_{k,L}),\mathrm{Im}(r^{[i]}_ky^{[i]}_{k,L})\Big]^T, \label{nextround}
\end{align}
where $y^{[i]}_{k,l}$ is the $l$-th element of $\mathbf{y}^{[i]}_k$.
It can be seen that $\mathbf{x}^{[i+1]}_k(0)$ is the starting point for the $(i+1)$-th iteration.
This completes one UMWFL round.

\section{Training Loss Minimization}

Ideally, the optimization of $\{\mathbf{F}^{[i]},r_k^{[i]},t_k^{[i]}\}$ should be performed to minimize the training loss, i.e., $\min~\Lambda(\mathbf{x}_k^{[i+1]}(0))$.
However, the analytical expression of $\mathbb{E}\left[\Lambda(\mathbf{x}_k^{[i+1]}(0))\right]$, where the expectation is taken over receiver noises and model parameters, is usually challenging to derive.
As a compromise approach, we resort the use of the upper bound of the expectation of $\mathbb{E}\left[\Lambda(\mathbf{x}_k^{[i+1]}(0))-\Lambda(\bm{\theta}^*)\right]$ as a metric to capture the degradation on training loss.
The upper bound depends on the MSE of the model parameters' estimation $\mathbb{E}\left[\|\mathbf{x}^{[i+1]}_k(0)-\bm{\theta}^{[i]}\|_2^2\right]$ where
\begin{align}
&\bm{\theta}^{[i]}=\sum_{k=1}^K\frac{|\mathcal{D}_k|}{\sum_{l=1}^{K}|\mathcal{D}_l|}\mathbf{x}_k^{[i]}(E) \label{global}
\end{align}
is equivalent to the gradient descent of the objective function of \eqref{FL}.
Specifically, based on \eqref{nextround} and \eqref{global}, the MSE between the received local parameter $\mathbf{x}^{[i+1]}_k(0)$ and the target local parameter $\bm{\theta}^{[i]}$ at the $i$-th UMWFL iteration is
\begin{align}\label{mse}
&\mathbb{MSE}_k^{[i]}\left(\mathbf{F}^{[i]}, r_k^{[i]},\{t_k^{[i]}\} \right)
=\mathbb{E}\left[\Big\|\mathbf{x}^{[i+1]}_k(0)-\bm{\theta}^{[i]}\Big\|_2^2\right]
\nonumber\\
&=2\eta^{[i]} S\Bigg[\gamma\sum_{j=1}^K\Big|r_k^{[i]}(\mathbf{g}_k^{[i]})^H\mathbf{F}^{[i]}\mathbf{h}_j^{[i]}\sqrt{t_j^{[i]}}-\alpha_j\Big|^2
+\gamma\sigma^2_b\|r_k^{[i]}(\mathbf{g}_k^{[i]})^H\mathbf{F}\|_2^2+\sigma_k^2|r_k^{[i]}|^2\Bigg],
\end{align}
where $\alpha_k=\frac{|\mathcal{D}_k|}{\sum_{l=1}^{K}|\mathcal{D}_l|}$ and the equality is due to \eqref{sk}, \eqref{Ri}, \eqref{yk}, \eqref{nextround}, \eqref{global}, and the independence among $\{\mathbf{s}_k|\forall k\}$.
Having obtained the MSE expression, the next step is to derive the relationship between the loss bound and the MSE.
To this end, we first introduce the following assumption imposed on the loss function.
\begin{assumption}
(i) The function $\Lambda(\mathbf{x})$ is $\mu$-strongly convex.
(ii) The function $\sum_{\mathbf{d}_{k,l}\in\mathcal{D}_k}\Theta(\mathbf{d}_{k,l}, \mathbf{x})$ is twice differentiable and satisfies $\sum_{\mathbf{d}_{k,l}\in\mathcal{D}_k}\nabla^2_{\mathbf{x}}\Theta(\mathbf{d}_{k,l}, \mathbf{x})\preceq L\mathbf{I}$.
\end{assumption}

Under Assumption 1, the relationship between $\Lambda(\mathbf{x}_k^{[i+1]}(0))$ and $\Lambda(\bm{\theta}^{*})$ is summarized in the following theorem.
\begin{theorem}
With $(\varepsilon,E)=(\frac{\sum_{k=1}^{K}|\mathcal{D}_k|}{KL}, 1)$, the UMWFL scheme satisfies
\begin{align}
&\mathbb{E}\left[\Lambda(\mathbf{x}_k^{[i+1]}(0))-\Lambda(\bm{\theta}^*)\right]
\leq
\sum_{i'=0}^{i}A^{[i']}\,\mathop{\mathrm{max}}_{k=1,\cdots,K}\mathbb{MSE}^{[i']}_k,
\end{align}
for any $\{\mathbf{F}^{[i']},r_k^{[i']},t_k^{[i']}\}_{i'=0}^i$ as $i\rightarrow+\infty$, where
\begin{align}
&A^{[i']}=
\frac{KL\left(3+2K^{-1}\right)}{2\sum_{k=1}^{K}|\mathcal{D}_k|}
\left(1-\frac{\mu\sum_{k=1}^{K}|\mathcal{D}_k|}{KL}\right)^{i-i'}.
\end{align}
\end{theorem}
\begin{proof}
We first derive the upper bound between consecutive parameter vectors.
Then the upper bound between consecutive losses can be derived.
Finally, we derive the upper bound of the sequence of loss function values.
Please refer to \cite{umair} for detailed derivations.
\end{proof}

Theorem 1 shows a diminishing $A^{[i']}\rightarrow 0$ for a large $i-i'$, meaning that the impact from earlier UMWFL iterations vanishes as the UMWFL continues.
On the other hand, if $\mathbb{MSE}^{[i']}_k\rightarrow 0$ for all $k$, then $\Lambda(\mathbf{x}_k^{[i+1]})$ is an unbiased estimate of $\Lambda(\bm{\theta}^*)$.
This demonstrates the effectiveness of UMWFL in the asymptotic regime.
The convexity and smoothness in Assumption 1 have been adopted in most loss bound analysis of FL (e.g., \cite{fed3}).
Although it seems to be restrictive for some realistic applications, analysis under Assumption 1 could provide important insights of the behavior of UMWFL in nonconvex cases.

Based on Theorem 1, the training loss minimization in UMWFL systems is formulated as
\begin{align}
\mathop{\mathrm{min}}_{\substack{\{\mathbf{F}^{[i']}\in\mathcal{F},r_k^{[i']},|t_k^{[i]}|^2\leq P_0\}}}
\quad&
\sum_{i'=0}^{i}A^{[i']}\,\mathop{\mathrm{max}}_{k=1,\cdots,K}\mathbb{MSE}^{[i']}_k.
\end{align}
It can be seen that the above problem and constraints can be decoupled for each iteration and the minimization at the $i$-th UMWFL iteration, $\forall i$, is given by
\begin{subequations}
\begin{align}
\mathcal{P}:\mathop{\mathrm{min}}_{\substack{\mathbf{F},\,\{r_k,t_k\}}}
\quad&
\mathop{\mathrm{max}}_{k=1,\cdots,K}~\gamma\sum_{j=1}^K\Big|r_k\mathbf{g}_k^H\mathbf{F}\mathbf{h}_jt_j-\alpha_j\Big|^2
+\gamma\sigma^2_b\|r_k(\mathbf{g}_k)^H\mathbf{F}\|_2^2+\sigma_k^2|r_k|^2
\label{P0}
\\
\quad\quad\quad\mathrm{s. t.}\quad\quad
&\mathbf{F}\in\mathcal{F},\quad |t_k|^2\leq P_0, \quad \forall k,
\end{align}
\end{subequations}
where the UMWFL iteration index $i'$ is removed since there is no dependence among different UMWFL iterations.

Problem $\mathcal{P}$ is generally NP-hard due to the unit-modulus constraints \cite{analog,analog2}.
In addition, the coupling between variables $\{r_k,t_k\}$ and $\mathbf{F}$ introduces nonlinearity and nonconvexity to the considered problem.
Finally, the large dimension of $\mathbf{F}$ call for the design of low-complexity algorithms in the scenario with massive numbers of antennas.

\section{Penalty Alternating Minimization for UMWFL}

In this section, the PAM algorithm, which consists of two nested layers of iterations (i.e., an outer-layer iteration and an inner-layer iteration), will be proposed to optimize the system performance.
Below we first introduce the outer-layer iteration.

\subsection{Outer-Layer Iteration Via Alternating Optimization}

To resolve the coupling between variables $\{r_k,t_k\}$ and $\mathbf{F}$, this paper adopts an alternating optimization framework \cite{beck}, which optimizes one design variable at a time with others being fixed.
Starting with an initial solution $\{\mathbf{F}^{(0)},r_k^{(0)},t_k^{(0)}\}$, the entire procedure
solving problem $\mathcal{P}$ for the $(n+1)$-th outer iteration, $\forall n$, can be elaborated below:
\begin{subequations}
\begin{align}
\mathbf{F}^{(n+1)}=\mathop{\mathrm{arg~min}}_{\substack{\mathbf{F}}}
&
\mathop{\mathrm{max}}_{k}~\Bigg(\sum_{j=1}^K\Big|r_k^{(n)}\mathbf{g}_k^H\mathbf{F}\mathbf{h}_jt_j^{(n)}-\alpha_j\Big|^2
+\sigma^2_b\|r_k^{(n)}\mathbf{g}_k^H\mathbf{F}\|_2^2\Bigg) \nonumber\\
\quad\quad\mathrm{s. t.}\quad
&|F_{l,l'}|=1,\quad \forall l,l', \label{problemF}
\\
\{r_k^{(n+1)}\}=\mathop{\mathrm{arg~min}}_{\substack{\{r_k\}}}
&
\mathop{\mathrm{max}}_{k}
\Bigg(
\gamma\sigma^2_b\|r_k\mathbf{g}_k^H\mathbf{F}^{(n+1)}\|_2^2+\sigma_k^2|r_k|^2
+
\gamma\sum_{j=1}^K\Big|r_k\mathbf{g}_k^H\mathbf{F}^{(n+1)}\mathbf{h}_jt_j^{(n)}-\alpha_j\Big|^2
\Bigg), \label{problemq}
\\
\{t_k^{(n+1)}\}=\mathop{\mathrm{arg~min}}_{\substack{\{t_k\}}}
&
\mathop{\mathrm{max}}_{k}~\sum_{j=1}^K\Big|r_k^{(n+1)}\mathbf{g}_k^H\mathbf{F}^{(n+1)}\mathbf{h}_jt_j-\alpha_j\Big|^2 \nonumber\\
\quad\quad\mathrm{s. t.}\quad
&|t_k|^2\leq P_0,\quad k=1,\cdots,K, \label{problemp}
\end{align}
\end{subequations}
where $\{\mathbf{F}^{(n)},r_k^{(n)},t_k^{(n)}\}$ is the solution at the $n$-th outer iteration.
The iterative procedure stops until $n$ reaches the maximum iteration number $n=N_{\mathrm{max}}$.

It can be seen that problem \eqref{problemq} is a typical least squares problem.
The optimal solution is found by setting the derivative $\partial \mathbb{MSE}_k/\partial \mathrm{conj}(r_k)$ to zero, which yields
\begin{align}
r_k^{(n+1)}=&
\frac{\sum_{j=1}^K\alpha_j\mathrm{conj}\left(\mathbf{g}_k^H\mathbf{F}^{(n+1)}\mathbf{h}_jt_j^{(n)}\right)}
{\sum_{j=1}^K|t_j^{(n)}\mathbf{g}_k^H\mathbf{F}^{(n+1)}\mathbf{h}_j|^2
+\sigma^2_b\|\mathbf{g}_k^H\mathbf{F}^{(n+1)}\|_2^2
+\frac{\sigma_k^2}{\gamma}}. \label{qk}
\end{align}
On the other hand, problem \eqref{problemp} is a convex optimization problem, which can be solved with a complexity of $\mathcal{O}(K^{3.5})$ via CVX, a Matlab software package for solving convex problems based on interior point method (IPM).
Consequently, the key challenge is to tackle the nonconvex problem \eqref{problemF}.

While problem \eqref{problemF} can be transformed into a convex problem via semidefinite relaxation (SDR), the computational complexity for solving the SDR problem is at least $\mathcal{O}\left(N^{7}\right)$ \cite{opt1}.
Since $N$ is usually much larger than $K$, this method is not desirable.
In the following, a new algorithm termed PAM, which decomposes \eqref{problemF} into smaller subproblems that are solved by closed-form updates, is proposed for achieving both excellent performance and significantly lower computational complexities.

\subsection{Inner-Layer Iteration Via PAM}

Since $\mathbf{F}$ is a matrix, its vectorization is given as $\mathbf{f}=\mathrm{vec}\left(\mathbf{F}\right)\in\mathbb{C}^{N^2\times 1}$.
Applying $\mathrm{Tr}\left(\mathbf{A}\mathbf{X}\mathbf{B}\mathbf{X}^T\right)=\mathrm{vec}(\mathbf{X})^T\left(\mathbf{B}^T\otimes\mathbf{A}\right)\mathrm{vec}(\mathbf{X})$ \cite{vec},
we have
\begin{align}
r_k^{(n)}\mathbf{g}_k^H\mathbf{F}\mathbf{h}_jt_j^{(n)}&
=(\mathbf{a}_{k,j}^{(n)})^H\mathbf{f},
\\
\sigma^2_b\|r_k^{(n)}\mathbf{g}_k^H\mathbf{F}\|_2^2&
=\mathbf{f}^H\mathbf{G}_k^{(n)}\mathbf{f},
\end{align}
where
\begin{align}
\mathbf{G}_k^{(n)}&=\sigma^2_b|r_k^{(n)}|^2\mathbf{I}_N\otimes\left(\mathbf{g}_k\mathbf{g}_k^H\right),
\\
\mathbf{a}_{k,j}^{(n)}&=\left[r_k^{(n)}t_j^{(n)}\left(\mathbf{h}_j^T\otimes\mathbf{g}_k^H\right)\right]^H.
\end{align}

Problem \eqref{problemF} is thus re-formulated as
\begin{align}
\mathcal{P}_F:\mathop{\mathrm{min}}_{\substack{\mathbf{f}}}
\quad&\mathop{\mathrm{max}}_{k=1,\cdots,K}~\left(
\sum_{j=1}^K\Big|(\mathbf{a}_{k,j}^{(n)})^H\mathbf{f}-\alpha_j\Big|^2+
\mathbf{f}^H\mathbf{G}_k^{(n)}\mathbf{f}\right)\nonumber\\
\quad\quad\quad\mathrm{s. t.}\quad
&|f_{l}|=1,\quad l=1,\cdots,N^2. \label{PF}
\end{align}
To handle the nonseparable objective function, variable splitting of $\mathbf{f}$ is proposed such that $\mathbf{f}=\mathbf{u}_1=\cdots=\mathbf{u}_K$, where $\{\mathbf{u}_k\}$ are auxilliary variables.
Moreover, to handle the unit-modulus constraints, another auxilliary variable $\mathbf{z}=\mathbf{f}$ is introduced.
For the newly introduced equality constraints, they can be transformed into quadratic penalties in the objective function \cite{penalty}.
Thus, $\mathcal{P}_F$ is approximately transformed into
\begin{align}
\mathop{\mathrm{min}}_{\substack{\mathbf{f},\mathbf{z},\{\mathbf{u}_{k}\}}}
\quad&\mathop{\mathrm{max}}_{k=1,\cdots,K}~\left(
\sum_{j=1}^K\Big|(\mathbf{a}_{k,j}^{(n)})^H\mathbf{u}_{k}-\alpha_j\Big|^2+
\mathbf{u}_{k}^H\mathbf{G}_k^{(n)}\mathbf{u}_{k}\right)
+\rho\left(\frac{1}{K}\sum_{j=1}^K\|\mathbf{u}_{j}-\mathbf{f}\|_2^2+\|\mathbf{z}-\mathbf{f}\|_2^2\right)
\nonumber\\
\mathrm{s. t.}\quad
&|z_{l}|=1,\quad l=1,\cdots,N^2, \label{PFa}
\end{align}
where $\rho$ is a penalty parameter for fine-tuning.
It can be proved that $\mathcal{P}_F$ and \eqref{PFa} are equivalent problems as $\rho\rightarrow+\infty$ \cite{penalty}.
However, this case also leads to the gradient norm of the objective function of \eqref{PFa} being infinite, making \eqref{PFa} difficult to solve.
Therefore, $\rho$ controls the tradeoff between approximation error and difficulty in solving \eqref{PFa}.

We address \eqref{PFa} using alternating minimization, in which the cost function is iteratively minimized with respect to one variable whereas the others are fixed.
Starting with an initial $\mathbf{f}^{(0)}=\mathbf{z}^{(0)}=\mathbf{u}_{k}^{(0)}=\mathrm{vec}\left(\mathbf{F}^{(n)}\right)$, the whole process consists of iteratively solving
\begin{subequations}
\begin{align}
\mathbf{u}_{k}^{(m+1)}=\mathop{\mathrm{arg~min}}_{\substack{\mathbf{u}_{k}}}
\quad&\sum_{j=1}^K\Big|(\mathbf{a}_{k,j}^{(n)})^H\mathbf{u}_{k}-\alpha_j\Big|^2+
\mathbf{u}_{k}^H\mathbf{G}_k^{(n)}\mathbf{u}_{k}
+\frac{\rho}{K}\|\mathbf{u}_{k}-\mathbf{f}^{(m)}\|_2^2,\quad \forall k, \label{PFa1}
\\
\mathbf{f}^{(m+1)}=\mathop{\mathrm{arg~min}}_{\substack{\mathbf{f}}}
\quad&\rho\Bigg(\frac{1}{K}\sum_{j=1}^K\|\mathbf{u}_{j}^{(m+1)}-\mathbf{f}\|_2^2
+\|\mathbf{z}^{(m)}-\mathbf{f}\|_2^2\Bigg), \label{PFa2}
\\
\mathbf{z}^{(m+1)}=\mathop{\mathrm{arg~min}}_{\substack{|z_{l}|=1, \forall l}}
\quad&\rho\|\mathbf{z}-\mathbf{f}^{(m+1)}\|_2^2, \label{PFa3}
\end{align}
\end{subequations}
where $m$ is the inner iteration index.
It can be verified that the objective function of \eqref{PFa} is strongly convex.
Therefore, despite the non-differentiability of the objective, the alternating minimization \eqref{PFa1}--\eqref{PFa3} is guaranteed to converge to a stationary point of \eqref{PFa} \cite{beck}.
The iterative procedure stops until $m$ reaches the maximum iteration number $m=M_{\mathrm{max}}$.

The remaining question is how to solve \eqref{PFa1}--\eqref{PFa3} optimally.
We notice that problems \eqref{PFa1} and \eqref{PFa2} are standard least squares problems, thus their solutions are given by the following closed-form expressions
\begin{align}
\mathbf{u}_{k}^{(m+1)}&=\left(
\sum_{j=1}^K\mathbf{a}_{k,j}^{(n)}(\mathbf{a}_{k,j}^{(n)})^H+\mathbf{G}_k^{(n)}+\rho\mathbf{I}
\right)^{-1}\left(\sum_{j=1}^K\alpha_j\mathbf{a}_{k,j}^{(n)}+\frac{\rho}{K}\mathbf{f}^{(m)}\right), \label{ukm}
\\
\mathbf{f}^{(m+1)}&=\frac{1}{2}\left(\frac{1}{K}\sum_{j=1}^K\mathbf{u}_{j}^{(m+1)}+\mathbf{z}^{(m)}\right), \label{fm}
\end{align}
respectively.
On the other hand, problem \eqref{PFa3} is the projection of $\mathbf{f}^{(m+1)}$ onto unit-modulus constraint and the optimal solution is simply
\begin{align}
\mathbf{z}^{(m+1)}=\mathrm{exp}\left(\mathrm{j}\angle\mathbf{f}^{(m+1)}\right). \label{zm}
\end{align}

\begin{figure*}[!t]
\centering
\subfigure[]{\includegraphics[width=0.99\textwidth]{map.eps}}
\subfigure[]{\includegraphics[height=0.42\textwidth]{loss.eps}}
\subfigure[]{\includegraphics[height=0.42\textwidth]{performance.eps}}
  \caption{Comparison between the proposed and benchmark schemes when $N=8$ with $K=3$: a) Multi-vehicle heterogeneous datasets in ``Town 05'' map; b) Training loss versus the number of UMWFL iterations; c) Average precision at $\textrm{IoU}=0.5$ on the testing datasets at different vehicles. The axis of each vehicle scales from average precision (AP) $0\%$ to $100\%$ with a step size of $10\%$.}
\end{figure*}

\subsection{Summary and Complexity Analysis of PAM}

In summary, the complete PAM algorithm for solving problem $\mathcal{P}$ is summarized in \textbf{Algorithm 1}.
The computational complexity (CP) is dominated by the steps at lines 4--6 in Algorithm 1, of which the CPs are given by $\mathcal{O}(KN^2)$, $\mathcal{O}(KN^2)$, and $\mathcal{O}(N^2)$, respectively.
Therefore, the total CP of PAM is $\mathcal{O}(N_{\rm{max}}M_{\rm{max}}KN^2)$

\begin{algorithm}[!t]
    \caption{UMWFL with PAM}
        \begin{algorithmic}[1]
            \State  Initialize $\mathbf{F}^{(0)},\,\{r_k^{(0)},t_k^{(0)}\}$ for all $k$.
            \State \textbf{For} $n=0,\cdots,N_{\rm{max}}-1$:
            \State \ \ \textbf{For} $m=0,\cdots,M_{\rm{max}}-1$:
            \State \ \  \ \  Update $\{\mathbf{u}_k^{(m+1)}\}_{k=1}^K$ using \eqref{ukm} with CP $\mathcal{O}(KN^2)$.
            \State \ \   \ \ Update $\mathbf{f}^{(m+1)}$ using \eqref{fm} with CP $\mathcal{O}(KN^2)$.
            \State \ \   \ \ Update $\mathbf{z}^{(m+1)}$ using \eqref{zm} with CP $\mathcal{O}(N^2)$.
            \State \ \ \textbf{End~For}
            \State \ \ Update $\mathbf{F}^{(n+1)}$ as $\mathbf{F}^{(n+1)}=\mathrm{mat}(\mathbf{z}^{(M_{\rm{max}})})$.
            \State \ \ Update $\{r_k^{(n+1)}\}_{k=1}^K$ using \eqref{qk} with CP $\mathcal{O}(KN^2)$.
            \State \ \ Update $\{t_k^{(n+1)}\}_{k=1}^K$ using \eqref{problemp} with CP $\mathcal{O}(K^{3.5})$.
            \State \textbf{End~For}
            \State Output $\mathbf{F}^{\diamond}=\mathbf{F}^{(N_{\rm{max}})}$, $r_k^\diamond=r_k^{(N_{\rm{max}})}$, $t_k^\diamond=t_k^{(N_{\rm{max}})}$.
        \end{algorithmic}
\end{algorithm}

\section{Results and Discussions}

This section presents simulation results to verify the performance of the proposed scheme.
We consider the object detection task in autonomous vehicle (AV) systems \cite{zijian,fedcav2}.
This task involves three main steps: 1) distributed dataset generation and storage at each vehicle; 2) collecting labels from nearby road-side infrastructures (RSIs); 3) federated learning within a vehicle platoon.
The three steps are executed sequentially and may not occur in the same place.
The case of $K=3$ is simulated and the number of total UMWFL iterations is set to $15$.
The average precision at intersection of union (IoU) equal to 0.5 is used for performance evaluation.
Besides the proposed PAM algorithm, we also simulate a baseline scheme, which sets $\mathbf{F}=\mathbf{I}_N$ and optimize $\{t_k,r_k\}$ using \eqref{problemq} and \eqref{problemp}.

For dataset generation and storage, it can be realized via a physical-world testbed; but this involves high implementation costs.
Car Learning to Act (CARLA) \cite{carla} is a widely-accepted unreal-engine platform that provides complex urban driving scenarios and high 3D rendering quality such that the AV object detection can be prototyped in virtual-reality.
Hence, in this paper, all the datasets are generated by CARLA.
In particular, we employ CARLA to generate $29$ vehicles in the ``Town05'' map, among which $3$ are autonomous driving vehicles that can generate the point-cloud data at a frequency of $10$ frames/s.
The entire dataset consists of $600$ frames at each vehicle, where $100$ frames are used for training and $500$ frames are used for testing.
Fig.~2a illustrates the bird eye view of the simulated world and the locations of all vehicles.
The datasets at different vehicles are non-i.i.d. due to various fields of views.
For label collection, it is assumed that each vehicle can obtain the ground truth label from its nearby RSI.
The labels should satisfy (Karlsruhe Institute of Technology and Toyota Technological Institute) KITTI formats.
This can be realized by processing the raw data generated from CARLA using the Python scripts in \cite{zijian,shanfeng}.
Finally, to simulate the vehicle platoon federated learning, we need to specify the local model and the communication parameters.
Specifically, the sparsely embedded convolutional detection (SECOND) neural network \cite{second} is adopted for object detection.
Each round of local training processes $100$ frames and the Adam optimizer is adopted with a learning rate of $10^{-4}$.
The number of local updates is $E=1$.
The model training is implemented by PyTorch using Python 3.8 on a Linux server with an NVIDIA RTX 3090 GPU.
For the communication model, the pathloss of the user $k$ is set to $\varrho_{k}=-40\,\mathrm{dB}$, and $\mathbf{h}_{k}$ and $\mathbf{g}_{k}$ are generated according to $\mathcal{CN}(\mathbf{0},\varrho_{k}\mathbf{I}_N)$.
The edge server is equipped with $N=8$ antennas.
The power scaling factor $\gamma=1$ and the maximum transmit powers at users are $P_0=1~\mathrm{W}$ (i.e., $30\,\mathrm{dBm}$).
The noise powers at the server and users are set as $-80\,\mathrm{dBm}$, which capture the effects of thermal noise, receiver noise, and interference.

\textbf{Performance}.
Firstly, it can be seen from Fig.~2b that the training loss of the baseline scheme fluctuates violently at the iterations $7$ and $11$, because the model errors are non-negligible due to the strong multi-path fading at these communication rounds.
In contrast, via proper phase shifting at the server, such a multi-path fading is effectively alleviated by the proposed scheme.
Therefore, the training loss of the proposed UMWFL with the PAM algorithm converges to a significantly smaller value than that of the baseline scheme.
Secondly, as shown in Fig.~2c, the average precision region of the proposed PAM scheme encloses that of the baseline scheme, which implies that the PAM scheme achieves higher learning performance on all the testing datasets.
Particularly, the performance improvement is between $10\%$ and $50\%$.
Lastly, the qualitative performance of the proposed and baseline schemes are compared in Fig.~3.
It can be observed that the proposed UMWFL with PAM detects all the objects correctly.
However, the baseline scheme misses a turning car at vehicle~$1$, misses an occluded car at vehicle~$2$, and outputs an inaccurate detection at vehicle~$3$.
This demonstrates the non-robustness of wireless FL if no phase optimization is performed.
In contrast, if non-unit-modulus beamforming optimization is employed, significantly higher implementation costs are involved.
The proposed UMWFL can effectively strike a balance between the robustness of the FL and the cost at the server, which is a competitive technology in real edge FL implementations.

\begin{figure*}[!t]
\centering
 \includegraphics[width=0.99\textwidth]{federated_all.eps}
  \caption{Detection results when $N=8$ with $K=3$. The red box is the ground truth; the blue box is from the proposed UMWFL scheme; the green box is from the fixed beamforming (baseline) scheme. Vehicle 1: The baseline scheme misses a turning car; Vehicle 2: The baseline scheme misses an occluded car; Vehicle 3: The baseline scheme outputs an inaccurate detection.}
\end{figure*}

\section{Conclusion}

This paper proposed the UMWFL scheme to support simultaneous transmission of local model parameters.
The training loss bound was derived and a low-complexity large-scale optimization algorithm was proposed to minimize the training loss.
The performance of the UMWFL framework with the proposed optimization algorithm was verified by using the CARLA autonomous driving platform.
It was found that the UMWFL scheme achieves robust FL performance and low communication costs.

\section{Acknowledgement}

This work was supported in part by the National Natural Science Foundation of China under Grant 62001203, in part by the Shenzhen Science and Technology Program under Grant RCB20200714114956153, and in part by the Shenzhen Fundamental Research Program under Grant JCYJ20190809142403596.


\section{Introduction}
Let $F$ be a finite field of characteristic $p$, where $p$ is a prime. An {\em additive code} of length $n$ over $F$ is a subset of $F^n$ with the property that for all $u,v \in C$, we have that $u+v \in C$. It is easy to prove that an additive code over $F$ is linear over ${\mathbb F}_p$, the finite field with $p$ elements. Thus, an additive code over $F$ can be defined as the row space over ${\mathbb F}_p$ of a matrix whose elements are from $F$. The code equivalence problem for additive codes is the following. Given two $k \times n$ matrices over $F$, when is there a permutation of the columns of one of the matrices, together with a permutation $\sigma_i$, for $i \in \{1,\ldots,n\}$, which is a permutation of the elements of $F$ in the $i$-th coordinate, so that the row-spaces over ${\mathbb F}_p$ are the same. If $F={\mathbb F}_p$ then an additive code is linear and if, furthermore, $p=2$ then the code is a binary linear code. If the permitations $\sigma_i$ are all additive then the codes are {\em additively equivalent}. It is not known if equivalent additive codes are necessarily additively equivalent, see \cite{BD2021}.

A multi-graph is a graph whose edges are assigned a weight from $\{1,\ldots,h\}$, for some natural number $h$. An edge joining vertices $u$ and $v$ of weight $w$ is interpreted as $w$ edges joining $u$ and $v$. The multi-graph isomorphism problem is the following. Given two multi-graphs determine if there a bijection between the set of vertices which induces a bijection on the edges. i.e. edges of weight $w$ are mapped to edges of weight $w$. If $h=1$ then a  multi-graph is simply a graph.

In \cite{PR1997}, Petrank and Roth provide a polynomial time reduction from the graph isomorphism problem to the binary linear code equivalence problem. The graph isomorphism problem is in NP but is not known to be NP-complete. The problem is not known to be solvable in polynomial time either and is therefore a good candidate to belong to the computational complexity class NP-intermediate.
Known algorithms for graph isomorphism include McKay's Nauty algorithm \cite{McKay1981}, Ullmann's algorithm \cite{Ullmann1976}, the VF2 algorithm \cite{CFSV2004} and the parameterised matching algorithm \cite{MKEHHP2013}. All these algorithms have exponential worst case performance. Solving isomorphism generally takes much longer time if there is no match, since all possible mappings are eventually searched until it is shown that there is no isomorphism. The latter extends to multi-graph isomorphism and is based on a parameterised sequence which is a walk that covers every vertex in the graph.

In this note we extend the polynomial time reduction of Petrank and Roth \cite{PR1997} to a polynomial time reduction of multi-graph isomorphism to additive code equivalence, where the code is over an extension of ${\mathbb F}_{2}$. Since a graph on $n$ vertices has at most $\frac{1}{2}n(n-1)$ edges, we assume that the multi-graph has $h \leqslant \frac{1}{2}n(n-1)$ weights and that these are from the set $\{1,\ldots, h\}$.

\section{The reduction of multi-graph isomorphism to additive code equivalence}

Let $\Gamma$ be a multi-graph with vertex set $V$, whose edges have weights belonging to the set $\{1,\ldots,h\}$. Let $A$ be the incidence matrix of $\Gamma$, whose rows are indexed by the edges, whose columns are indexed by the vertices and where the edge-vertex entry is equal to the weight of the edge if the edge is incident with the vertex and zero otherwise. Let $E$ denote the set of edges of $\Gamma$ and let $E_i$ be the subset of $E$ of edges of weight at least $i$ (so $E_1=E$).

Let $e$ be a primitive element of ${\mathbb F}_{2^r}$, where $r>\log_2 h$. Let $\mathrm{D}_i$ be the $|E_i| \times |E_i|$ diagonal matrix indexed by the edges of $E_i$ whose diagonal entry is $e^{j-1}$, where $j$ is the weight of the edge indexing the row of $\mathrm{D}_i$.

Let
$$
N=(h+2)|E_1|+|E_2|+\cdots+|E_{h}|+|V|.
$$
We map $\Gamma$ to the additive code which is the ${\mathbb F}_2$-row span of the $|E| \times N$ matrix
$$
\mathrm{G}(\Gamma)=\left( \underbrace{D_1 \ | \ D_1 \ | \ \cdots \ | \ D_1}_{h+2} \ | \ \begin{array}{c} O_2 \\ D_2 \end{array}
\ | \ \begin{array}{c} O_3 \\ D_3  \end{array} \ | \ \cdots \ | \ \begin{array}{c} O_{h} \\ D_{h}  \end{array}
\ | \ A
\right)
$$
where $O_i$ is a matrix of zeros, whose dimensions are determined by the fact that the matrix $\mathrm{G}(\Gamma)$ is a $|E| \times N$ matrix and the matrix $\mathrm{D}_i$ is a $|E_i| \times |E_i|$ matrix.

The following two observations will be vital.
\begin{enumerate}
\item[(O1)]
The row of $G(\Gamma)$ indexed by an edge of weight $i$ is a codeword of weight $h+3+i$ whose coordinates are either $0$ or $e^{i-1}$.
\item[(O2)]
The only codewords of weight at most $2h+3$ are rows of $G(\Gamma)$.
\end{enumerate}

The following is the main theorem of this note.

\begin{theorem} \label{mainthm}
The multi-graphs $\Gamma$ and $\Gamma'$ are isomorphic if and only if the additive codes generated by $\mathrm{G}(\Gamma)$ and $\mathrm{G}(\Gamma')$ are equivalent.
\end{theorem}

\begin{proof}

Suppose the multi-graphs $\Gamma$ and $\Gamma'$ are isomorphic. Then there is a permutation of the columns of $\mathrm{A}$ and a permutation of the rows of $\mathrm{A}$, which gives $A'$, the incidence matrix of the multi-graph $\Gamma'$. Apply the column permutation to the last $|V|$ columns of $\Gamma$ and the permutation of the rows of $\mathrm{A}$ to the rows of $\mathrm{G}(\Gamma)$. Since the permutation of the rows of $\mathrm{A}$, takes edges of weight $i$ to edges of weight $i$, it takes codewords of weight $h+3+i$, which are rows of $G(\Gamma)$, to codewords of weight $h+3+i$, which are rows of $G(\Gamma')$ by (O1) and (O2). Thus, $G(\Gamma)$ and $G(\Gamma')$ generate equivalent codes.

Now suppose the additive codes generated by $\mathrm{G}(\Gamma)$ and $\mathrm{G}(\Gamma')$ are equivalent. Then, by definition, there is an $N \times N$ permutation matrix $P$, a $|E| \times |E|$ change of basis matrix $S$ and permutations $\sigma_i$, $i\in \{1,\ldots,N\}$, such that we can apply $P$ to the columns of $\mathrm{G}(\Gamma)$, $S$ to the rows of $\mathrm{G}(\Gamma)$ and the permutation $\sigma_i$ to the $i$-th coordinate and obtain a matrix whose generates, over ${\mathbb F}_2$, the same code as $\mathrm{G}(\Gamma')$.

Properties (O1) and (O2) imply that the change of basis matrix $S$ is in fact a permutation matrix, since rows of $\mathrm{G}(\Gamma)$ must be mapped to rows of $\mathrm{G}(\Gamma')$.
Furthermore, $S$ maps codewords of weight $h+3+i$ to codewords of weight $h+3+i$, so the initial $N-|V|$ columns of $SG$ contain each vector of weight one, with a non-zero coordinate indexed by an edge of weight $i$, precisely $h+1+i$ times. Hence, the first $N-|V|$ columns of $G'$ can be obtained by permuting the first $N-|V|$ columns of $SG$. This permutation of the columns then maps the columns of $A$ to columns of $A'$.

The permutation $\sigma_i$ must fix all elements of $F$ appearing in the $i$-th column of a column of $A$, since by property (O1) the rows of $\mathrm{G}(\Gamma')$ of weight $h+3+i$ have coordinates $0$ or $e^{i-1}$.

Hence, we obtain $\mathrm{G}(\Gamma')$ from $\mathrm{G}(\Gamma)$ by applying a permutation of the columns of $A$ and a permutation of the rows of $A$. Thus, $\Gamma$ and $\Gamma'$ are isomorphic.
 \end{proof}

The polynomial time reduction follows from Theorem~\ref{mainthm}, by noting that $|N|<chn^2<n^4$, for some constant $c$.
Note also that we can replace equivalence with additive equivalence, since the permutations $\sigma_i$ are always trivial.

It is interesting to ask if $N$ could be decreased, possibly by increasing the size of the field extension. For $h=2$ it is possible to reduce $N$ to $3|E_1|+|E_2|+|V|$ by simply deleting one of the initial $D_1$ matrices. One can check that the proof still works, although a little more subtlety is required. Observe that the sum of two codewords of weight $5$ (corresponding to edges of weight one) cannot produce a codeword of weight $6$ since this would require a repeated single edge, which is an edge of weight two.

\section{Code equivalence algorithms}

We discussed known algorithms for solving the graph isomorphism problem in the introduction. Algorithms for solving binary linear code equivalence include Bouyukliev's algorithm \cite{Bouyukliev2007}, which is similar to McKay's graph isomorphism algorithm \cite{McKay1981}, Leon's algorithm \cite{Leon1982} and the support splitting algorithm of Sendrier \cite{Sendrier1999}, \cite{Sendrier2000}. The latter does not apply to all binary linear codes but it is interesting because it is fast for binary linear codes in which the dimension of the hull $H(C)=C \cap C^{\perp}$ is small. It seems that this is the case for codes obtained from the polynomial time reduction of the graph isomorphism problem due to Petrank and Roth \cite{PR1997}. This warrants further investigation.

The support splitting algorithm allows one to determine signatures from the weight distribution of the hull of truncated binary linear codes, truncating one coordinate at a time.  This works because one of the hulls of either the truncated code or the truncations of its dual is non-trivial. This can be seen as follows. Let $C_i$ be the code obtained from $C$ be setting the $i$-th coordinate to zero in all codewords. Suppose that $H(C)=\{0\}$. Thus, ${\mathbb F}_2^n=C\oplus C^{\perp}$.

Without loss of generality, we consider a truncation on the first coordinate. We have that
$$
(1,0,\ldots,0)=u+v,
$$
for some $u \in C$ and $v \in C^{\perp}$.

If $(u_1,v_1)=(1,0)$ then, since $v_1=0$, we have that $v \in (C_1)^{\perp}$ and since $v= (1,0,\ldots,0)+u$, we have that $v \in C_1$. Hence, $v \in H(C_1)$.

Similarly, if $(u_1,v_1)=(0,1)$ then, since $u_1=0$, we have that $u \in ((C^{\perp})_1)^{\perp}$ and since $u= (1,0,\ldots,0)+v$, we have that $u \in (C^{\perp})_1$. Hence, $u \in H((C^{\perp})_1)$.

This truncation trick does not carry over to additive codes. It is possible that both $H((C^{\perp})_1)$ and $H(C_1)$ are trivial. However, if we are interested in establishing additive equivalence then we can employ a slightly modified support splitting algorithm; instead of removing the $i$-th coordinate we take a subspace. In effect, truncating in the binary case would be equivalent to taking the subspace $\{0\}$. This is done in the following way. The $i$-th column of a generator matrix for $C$ is a vector $v_i \in {\mathbb F}_{2^r}^k$, where $|C|=2^k$. Writing $v_i$ out over the basis $\{1,\alpha,\ldots,\alpha^{r-1}\}$ of ${\mathbb F}_{2^r}$ over ${\mathbb F}_2$,
$$
v_i=\sum_{j=0}^{r-1} v_i^{(j)}\alpha^j,
$$
where $v_i^{(j)} \in {\mathbb F}_{2}^k$. Let $\rho_i$ be the subspace of  ${\mathbb F}_{2}^k$, of dimension at most $r$, spanned by $v_i^{(0)},\ldots, v_i^{(r-1)}$. The subspace $\rho_i$ is unaffected by an additive permutation on the $i$-th coordinate, which only has the effect of changing the basis. Thus, we can select a subspace $\pi_i$ of $\rho_i$ and replace the $i$-th column of the generator matrix with the corresponding vector of ${\mathbb F}_{2^r}^k$ obtained from the subspace $\pi_i$ by reversing the construction of $\rho_i$ above. This code we denote by $C_{i,\pi_i}$, where $i$ is the selected coordinate and $\pi_i$ is the subspace. For equivalent codes, we have that selecting a coordinate and a subspace in $C$, there must be a corresponding coordinate and subspace in $C'$ which is equivalent, i.e. $C_{i,\pi_i}$ is equivalent to $C'_{i',\pi'_{i'}}$, for some $i'$ and $\pi'_{i'}$. The weight distribution of the hull of $C_{i,\pi_i}$ can then be used as a signature with an aim of establishing the permutation taking $C$ to $C'$, as in the support splitting algorithm.