\section{Introduction}
\footnotetext[1]{Our code is available on Github at \url{https://github.com/Derek-Wds/MAD-VAE}}%
With the development of deep learning in recent years, research about image classification have made significant breakthroughs. Models such as ResNet \cite{xie2017aggregated} and VGG \cite{simonyan2014very} as well as their variants have achieved amazing accuracy on image classification tasks. Despite their high accuracy on image classification, neural network based models are often vulnerable to misclassification given a well-crafted input by an adversary. The differences between these \textit{adversarial examples} and original data are often imperceptible to humans, but result in misclassification by the model \cite{goodfellow2014explaining, szegedy2013intriguing}.

Recent research has considered two main attack methods: White-Box attack and Black-Box attack. Under the White-Box attack scenario, the attackers have access to, at the very least, the classification model as well as the pretrained parameters. Under the Black-Box attack situation, the attackers do not have any information about the classification model \cite{xu2019adversarial}.

Along with attacks, various defenses mechanisms have also been proposed in recent years. The defenses methods can mainly be divided into three groups \cite{xu2019adversarial}: (1) obfuscating the gradient of the classification model \cite{defense-vae18, samangouei2018defense, goodfellow2014explaining, kurakin2016adversarial, dong2017boosting}, (2) robust optimization such as such as training the classification model with adversarial examples \cite{goodfellow2014explaining} and (3) adversary detection before feeding images to the classification model \cite{grosse2017statistical, Xu_2018}.

In our research, we focus on the first kind of defense mechanism, gradient obfuscation, and improve upon the Defense-VAE \cite{defense-vae18} model. The idea of Defense-VAE is to use VAE's ability to learn the distribution of input data in order to transfer the potentially adversarial input data back to the data manifold. In the original VAE framework \cite{kingma2013auto}, the model consists of two parts: an encoder $E$ which maps the input data onto a latent distribution, and a decoder $D$ which tries to reconstruct the input data based on the variable $z$ sampled from the latent distribution. Defense-VAE leverages this idea in their defense model and proves that it is a robust way to prevent against attacks. Unfortunately, although generative defense models such as Defense-GAN \cite{samangouei2018defense} and Defense-VAE \cite{defense-vae18} have achieved decent results on prevent certain attacks, they have been shown to still be vulnerable to carefully crafted adversarial examples that attack the latent space of the generative model instead of the classifier directly \cite{jalal2017robust}.

The following paper consists of following parts: for theoretical background and related works, please see Appendix \ref{related}; in Section \ref{proposed}, we will introduce the research work we are focusing on; Section \ref{experiments} will be about the experiments and results while Section \ref{conclusion} discusses future works and concludes the paper. There are additional plots and tables in the remainder of the Appendix that were not included in the main body of the paper.

\section{Proposed MAD-VAE} \label{proposed}
We propose modifying the training procedure and introduce several additional loss functions, which help Defense-VAE to be more robust against adversarial attack \cite{jalal2017robust}. We call our model Manifold Awareness Defense-VAE because our new training steps attempt to capture the topological structure of the data manifold as well as the latent variables. In the following subsections, we will first talk about the motivation of our research and then describe our methods.

\subsection{Motivation}
In order to further mitigate the problem of adversarial attack, we formulate our objective in the following way: we want the encoded $z$ to represent the input data well while having the decoded $G(z)$ be correctly classified. Since $z_i$ follows a certain distribution for each class $i$, we want the points $z_i^j$ sampled from this distribution to be close to each other while maintaining that this cluster is well separated from clusters of other classes. Ideally, our model should be able to encode the adversarial data $X'$, where the benign $X$ data has label $i$, onto or close to the distribution $z_i$. Therefore, when we sample from the latent distribution and try to find $z_i^j$ that decodes to similar data, not only will it be less likely to decode data that lies in the other clusters but also the correct classification can be ensured.

\subsection{Specification}
We introduce three loss functions that increase the robustness of Defense-VAE through improving the awareness of the underlying topology of the data manifold. This is a realization of the goals mentioned in the previous section.

\subsubsection{Classification Loss}
When the classifier gets an input data $X$, it will extract certainly features from $X$ and make classification based on them. Since our goal is to make sure that the output of our model can be correctly classified, one straight forward way to improve the model is to having a classification loss on the output of the VAE $G(z)$. This idea is from GAN (Generative Adversarial Networks) \cite{goodfellow2014generative} where we have a generative model $G$ to generate data and a discriminator $D$ to make sure the generated data is what we want. A similar idea was proposed by Jang et al. \cite{jang2019need} for improving the topological awareness of generative defense models. Specifically, in our case, we have the VAE as our generative model and the pretrained classifier to be our discriminator model.

During training of the VAE, we fix the parameters of the classifier and feed the output of VAE $G(z)$ to the classifier to get a classification loss $L_c$. The new loss function of our model is as follows:
\begin{align*}
    L = L_r + \beta \cdot L_{kl} + \alpha \cdot L_c
\end{align*}
where $L_r$ is the reconstruction loss and $L_{kl}$ is the KL-divergence loss in the \textit{ELBO} equation. $\alpha$ and $\beta$ are two hyperparameters of the weight on different loss.

By doing so, we hypothesize that the classification loss can help to coerce the VAE model to generate features that are sensitive to the classifier, and at the same time, force the latent $z$ to lie in the cluster of $z_i$ for each class $i$.

\subsubsection{Proximity and Distance Loss}

According to Zhao et al. \cite{zhao2017deeper}, due to the nature of MNIST, it would be hard for a model like our own to learn a clear separation between classes and as a result the latent distributions of $z$ will certainly mix with each other. In order to solve this problem, we decide to add another explicit loss to constrain the latent variables. Similar to the idea from the research conducted by Mustafa et al. \cite{Mustafa_2019_ICCV}, we try to learn the center $c_i$ of the cluster of $z_i$ for each class $i$ and, at the same time, keep the clusters away from each other. The proposed loss of function can be described as follows:
\begin{align*}
    L_p & = ||z_i - c_i||_2 \\
    L_d & = \frac{1}{k-1} \cdot \sum_{j\neq i} \left(||z_i - c_j||_2 + ||c_i - c_j||_2\right) \\
    L_{pd} & = \sum_i\left[\alpha \cdot L_p - \sigma \cdot L_d \right]
\end{align*}
where $k$ is total number of labels in the dataset and the distances are computed by the euclidean distances. $\alpha$ and $\sigma$ are the weights between the proximity loss and distance loss respectively. Therefore, the loss function for our VAE model becomes:
\begin{align*}
    L = L_r + \beta \cdot L_{kl} + L_{pd}
\end{align*}
where the variables are the same as those above. The idea of this loss function is to force the latent variable $z_i$ with the same label $i$ to gather together towards the center $c_i$ while make sure $c_i$ and $z_i$ is farther away from the center of cluster of $z_j$ where $j \neq i$. Through the adversarial training on our defensive model, we are pushing the decision boundaries between each cluster to better separate the classes.

\subsubsection{Combined Loss}
Since we have discussed two new loss functions for generative defense models in the subsections above, it would be natural to think of a way to combine them together in order to utilize both their advantages. A naive approach to the aggregation of the two defense loss functions would be to design a weighted sum. To do so, we use the combined loss below for our VAE model:
\begin{align*}
    L = L_r + \beta \cdot L_{kl} + \gamma \cdot L_c + L_{pd}
\end{align*}
where $L_r$ is the reconstruction loss and $L_{kl}$ is the KL-divergence loss in the $ELBO$ equation. $\beta$ and $\gamma$ are two hyperparameters controlling the weight assigned to the different loss functions.

\section{Experiments} \label{experiments}
We compare our methods against the original Defense-VAE model, however, it may be noted that our model does not use any of the Defense-VAE's code base. Through our reproduction of Defense-VAE, we were able to match the performance detailed in the Defense-VAE paper \cite{defense-vae18}, therefore validating its functional equivalence.

We evaluate our defensive model under the FGSM, Rand-FGSM and CW white-box attacks as well as three other attack methods not included in our training process: PGD \cite{aleks2017deep}, momentum iterative FGSM \cite{dong2017boosting} and single pixel attack \cite{narodytska2016simple}. We also evaluate MAD-VAE and Defense-VAE under the Overpowered Attack \cite{jalal2017robust}. Although the Overpowered Attack is technically a white-box attack, it is a \textit{latent space} attack and therefore it is not reasonable to compare it against other white-box attacks. Instead, the Overpowered Attack gets its own section, after Black-Box Attacks.

Our classifier and MAD-VAE implementation are written using the PyTorch \cite{paszke2019pytorch} library and we use the open source package AdverTorch \cite{ding2019advertorch} for performing various attack methods. The Overpowered Attack is also written in PyTorch and is largely based off of the code provided by Jalal et al. \cite{jalal2017robust}.

We use the MNIST dataset \cite{lecun-mnisthandwrittendigit-2010} for our experiments, which contains 60,000 training images and 10,000 testing images. We generate a large, combined training dataset using three attack methods: FGSM, Rand-FGSM and CW with three different parameter settings for each one of them. We then use another parameter setting to generate the validation dataset to evaluate the performance of our models, and we use the testing data set to get our test dataset for all the six attack algorithms mentioned above. Since the attack generating process is stochastic, it is rare for non-empty intersections of different dataset to occur.

Due to time limitations and the computational complexity of the models, we have not been able to fully reproduce the experiments conducted in Defense-VAE \cite{defense-vae18} and Defense-GAN \cite{samangouei2018defense}. Our defense model, however, remains based on the architecture of Defense-VAE and is able to match the Defense-VAE \cite{defense-vae18} accuracy without any of our additional loss functions. The structure of MAD-VAE can be found in Table \ref{table:defense-vae} in the Appendix. In our experiments, we choose the classifier MagNet \cite{meng2017magnet} to be our classification model. MagNet's model structure can be found in Table \ref{table:classifier} in the Appendix.

The original Defense-VAE model converges after $5$ epochs while our model, with the new loss functions the proximity and distance loss variant converges after $5$ epochs, while the other variants converge after $10$ epochs. In order to achieve the best performance, we validated our model on a validation data set for different parameter settings. From our experiments on the validation set, we chose the weight of our loss functions as: $\alpha=0.01$, $\sigma=0.00001$, $\beta=0.1$ and $\gamma=0.1$. The other parameters setting are kept the same as Defense-VAE paper. We use the Adam and SGD optimizers for our gradient descent with a exponential learning rate scheduler.

\subsection{Results on White-Box Attack}
In this section, we will present the results on White-Box attack using six different attack methods: FGSM, Rand-FGSM, CW, MI-FGSM, PGD and Single Pixel. All of the adversarial images are generated using the AdverTorch \cite{ding2019advertorch} package with default parameters. The results are shown in Table \ref{table:whitebox-result}.

\input{whitebox_table}

We note that the attacks are successfully fooling the classifier except for the Single Pixel attack. Table \ref{table:whitebox-result} demonstrated that our methods make the original Defense-VAE more robust under all the six attacks. Additionally the model utilizing Classification loss outperformed the others on most of the situations with the Combined loss outperforming under the other attacks. Unlike what was does in the Defense-VAE paper, we did not fine-tune the classifier based on the output of out model. Through our topologically aware training, we consistently outperform Defense-VAE on the tested white-box attacks.

\subsection{Results on Black-Box Attack}
In this section, we present the results of our experiments of MAD-VAE compared with the vanilla Defense-VAE on the FGSM Black-Box attack. As what is mentioned above, attackers will have no access to the classifier information under the Black-Box attack scenario. They thus training a substitute model for generating attacks. The detail of substitute models can be found in  Table \ref{table:black-box}.

In Table \ref{table:blackbox-result}, we present our classification results for different substitute classification model from the Defense-VAE paper \cite{defense-vae18}. For all the experiments, we set the FGSM parameter $\epsilon$ to be $0.3$ which is the default value given by the original paper. From the table, we can spot that the Black-Box FGSM attack can reduce the classification accuracy of different classifiers to only $5$ percent. All models, including Defense-VAE, perform well on the Black-Box FGSM attacks. However, the Classification loss model and Proximity \& Distance model outperform Defense-VAE by several percentage points.

\input{blackbox_table}

\subsection{Results under the Overpowered Attack}
The Overpowered Attack \cite{jalal2017robust} is a powerful adversarial attack that targets the latent space of generative models. In the paper by Jalal et al., the authors use the Overpowered Attack to reduce the accuracy of Defense-GAN to 3\%. Defense-VAE is not only supposed to be a faster alternative to the online-optimization of Defense-GAN, but also provides a more robust defense generative model. Because there is some noise induced by the generative defense models, attacks with different noise robustness parameters were performed both in Jalal et al. experiments and our experiments. Because this attack is a latent space attack, the classifier cannot be directly targeted to yield a baseline accuracy, hence we trained a standard  ``Identity'' VAE on MNIST.

We measured the performance of the defensive models under two metrics. In the paper by Jalal et al., they measured the accuracy of the classifier to be the number of digits such that the classifier correctly classified the digit under all noise robustness parameters, referred to as ``Original Accuracy'' in our results. On the other hand, we also measured the accuracy of the classifier on the combined pool of adversarial images, regardless of the noise robustness parameters, which we refer to as ``Accuracy''. The results from the Overpowered Attack experiments can be found in Table \ref{table:op-result}

\input{OP_table}

From Table \ref{table:op-result} we can see that under both our accuracy and the original accuracy metrics, the vanilla Defense-VAE outperformed all MAD-VAE variants, although the Proximity \& Distance variant came within 1\% of the accuracy of the vanilla Defense-VAE. These results are initially very startling due to the decrease in adversarial robustness of the MAD-VAE models in comparison to the Defense-VAE model – we expected that the MAD-VAE accuracy would either increase or remain the same relative to Defense-VAE's accuracy. The intuition for this is seeded in Section \ref{sec:clustering-op} in which we plot the distribution of these adversarial data points, both in latent and ambient space.

\subsection{Clustering for White-Box Attacks}
In this section, we will show the clustering of the adversarial examples, latent variables, and VAE output images using the Uniform Manifold Approximation and Projection (UMAP) \cite{2018arXivUMAP} dimensional reduction technique. Clustering of VAE output images can give us a direct sense of the classification result, where images with same label tend to be close to each other. Similarly, since we are trying to cluster variables in latent space, it is helpful to use clustering figures to help us evaluate the results. Since we have four models in total, in the following plots we will denote vanilla Defense-VAE as $A$, Defense-VAE with combined loss as $B$, Defense-VAE with proximity and distance loss as $C$, and Defense-VAE with classification loss as $D$.

Since the experiments in the previous section show the strong performance of our MAD-VAE model endowed with the classification loss functions, due to the limited space we decided to only show the clustering plot of that model in this section. Additional plots for other defensive models and adversarial attacks can be found in the appendix.

\begin{figure}[H]
    \centering
    \includegraphics[width=.45\linewidth]{images/MNIST.png}
    \caption{Clustering plots of the benign test dataset.}
    \label{fig:benign-data}
\end{figure}

Figure \ref{fig:benign-data} is the clustering of original MNIST test dataset, from which we are able to see that images of digit $6$ and $0$ are well separated from others although there have some points have similar characteristics and lie in the wrong clusters. The images of digit $9$ and $4$ are mixed with each other due to the similarities of certain features in the images, while other clusters are close to each other but still well separated.

\begin{figure}[H]
\begin{subfigure}{.33\textwidth}
  \centering
  \includegraphics[width=\linewidth]{images/fgsm_pure.png}
  \caption{FGSM attack.}
  \label{fig:fgsm-clustering}
\end{subfigure}%
\begin{subfigure}{.33\textwidth}
  \centering
  \includegraphics[width=\linewidth]{images/r-fgsm_pure.png}
  \caption{Rand-FGSM attacks}
  \label{fig:rfgsm-clustering}
\end{subfigure}
\begin{subfigure}{.33\textwidth}
  \centering
  \includegraphics[width=\linewidth]{images/cw_pure.png}
  \caption{CW attack.}
  \label{fig:cw-clustering}
\end{subfigure}
\caption{Clustering plots for adversarial examples}
\label{fig:clustering}
\end{figure}

Figure \ref{fig:clustering} is the MNIST clustering of data and predicted labels generated from three attacks: FGSM, Rand-FGSM and CW. From the figure, we can still see clustering of the data, however, the clustering does not correspond to classes like in \ref{fig:benign-data}.

In the following two groups of plots, we will show the clustering of the output of our Defense-VAE model with the classification loss function as well as the latent variables $z$. Since our goal is to achieve data manifold awareness during training, we are hoping to spot clustering in our data and latent variables.

Figure \ref{fig:output-clustering} shows the clustering result of the output from Defense-VAE with classification loss on adversarial data. The colored labels of Figure \ref{fig:output-clustering} are predicted using the classifier. From this figure, we are able to see that the output of our model is clearly clustered by color, whereas the adversarial input (Figure \ref{fig:clustering}) is not. Although some of the data with different labels are still mixed, it overall demonstrate the potential ability of the model to remap data back to the original data manifold.

\begin{figure}[H]
\begin{subfigure}{.33\textwidth}
  \centering
  \includegraphics[width=1.\linewidth]{images/fgsm_classification.png}
  \caption{FGSM attack.}
  \label{fig:fgsm-output-clustering}
\end{subfigure}%
\begin{subfigure}{.33\textwidth}
  \centering
  \includegraphics[width=1.\linewidth]{images/r-fgsm_classification.png}
  \caption{Rand-FGSM attacks}
  \label{fig:rfgsm-output-clustering}
\end{subfigure}
\begin{subfigure}{.33\textwidth}
  \centering
  \includegraphics[width=1.\linewidth]{images/cw_classification.png}
  \caption{CW attack.}
  \label{fig:cw-output-clustering}
\end{subfigure}
\caption{Clustering plots of output from MAD-VAE with classification loss across different attack methods.}
\label{fig:output-clustering}
\end{figure}

Figure \ref{fig:latent-clustering} shows the clustering of the latent variable $z$ in the classification MAD-VAE. Since we are adding the classification loss with respect to the latent variables in our new loss function, we would hope to see the clear separation of latent variables compared to the vanilla Defense-VAE.

\begin{figure}[H]
\begin{subfigure}{.33\textwidth}
  \centering
  \includegraphics[width=\linewidth]{images/latent_fgsm_classification.png}
  \caption{FGSM attack.}
  \label{fig:fgsm-latent-clustering}
\end{subfigure}%
\begin{subfigure}{.33\textwidth}
  \centering
  \includegraphics[width=\linewidth]{images/latent_r-fgsm_classification.png}
  \caption{Rand-FGSM attacks}
  \label{fig:rfgsm-latent-clustering}
\end{subfigure}
\begin{subfigure}{.33\textwidth}
  \centering
  \includegraphics[width=\linewidth]{images/latent_cw_classification.png}
  \caption{CW attack.}
  \label{fig:cw-latent-clustering}
\end{subfigure}
\caption{Clustering plots for latent variables of the MAD-VAE with classification loss given adversarial examples as input.}
\label{fig:latent-clustering}
\end{figure}

\subsection{Clustering for Overpowered Attack} \label{sec:clustering-op}
For the clustering of Overpowered Attacks, we feel like it would be most representative of the attack to plot only one noise robustness set at a time. Therefore, in this section we only present the clustering for the Overpowered Attack with a noise robustness term of $1.0$. Additionally, our plots will focus on the Overpowered Attack against the proximity and distance loss MAD-VAE, since the effect of the adversary is most pronounced under this model.

Figure \ref{fig:mnist-op} shows the clustering of the output generated by the Overpowered Attack when attacking MAD-VAE with proximity and distance loss. In contrast to Figure \ref{fig:output-clustering}, we see that there is no obvious clustering structure associated with the underlying MNIST dataset. This is due to the fact that instead of attacking the ambient space of the classifier, the Overpowered Attack instead attacks the latent space of the defensive model. Intuitively, if the latent space accurately models the underlying data manifold of MNIST, then latent variables that lie outside of the learned manifold do not necessarily have to follow the underlying clustering when decoded.

\begin{figure}[H]
    \centering
    \includegraphics[width=.45\linewidth]{images_op/mnist_proxi.png}
    \caption{Adversarial clustering under the Overpowered Attack on MAD-VAE with proximity and distance loss.}
    \label{fig:mnist-op}
\end{figure}

In the next two plots, we will look at the effect of the Overpowered Attack on the latent variables of the defensive models. For clarity, we have split the original data and the adversarial data plots into two separate plots, though they should be viewed together. Figure \ref{fig:latent-op-clustering} shows the latent embedding of both the original data and the adversarial data generated by the Overpowered attack.

\begin{figure}[H]
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.75\linewidth]{images_op/latent_proxi_orig.png}
  \caption{Unperturbed data clustering}
  \label{fig:op-orig}
\end{subfigure}%
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.75\linewidth]{images_op/latent_proxi.png}
  \caption{Overpowered Attack clustering}
  \label{fig:op-latent}
\end{subfigure}
\caption{Clustering plots for latent variables of the MAD-VAE with proximity and distance loss. First figure is original data and second under the Overpowered Attack.}
\label{fig:latent-op-clustering}
\end{figure}

What we can see from Figure \ref{fig:latent-op-clustering} is that the Overpowered Attack tries to pull the original data points towards other clusters of data. The decoded images of the perturbed data are therefore more similar to other classes of numbers. The moving out data out of well-defined clusters towards other points could explain why Defense-VAE outperforms MAD-VAE in these tests, as Defense-VAE has a less structured latent space, as is shown in the Appendix Figure \ref{fig:vanilla-latent-clustering}. This, however, remains pure conjecture.

\subsection{Comments on the Overpowered Attack}
By exploiting the latent space of defensive models, the Overpowered Attack \cite{jalal2017robust} is very effective at finding what traditional literature on the subject considers adversarial examples. The traditional methodology surrounding adversarial examples is to find some adversarial perturbation under some $\ell_2$ distance budget from the original data point. It is natural to limit adversarial examples in this way since it prevents examples from becoming too noisy and defeating the ``imperceptible'' nature of adversarial examples.

While the methodology of limiting perturbations under some $\ell_2$ distance makes sense for ambient space attacks, for latent space attacks on defensive models, this methodology falls apart. The Overpowered Attack is one adversarial example generation method that seems to abuse this $\ell_2$ constraint. How does one abuse a distance metric? Well, latent space attacks abuse how inadequately the $\ell_2$ norm represents perceptual distance. The problem imposed by the $\ell_2$ norm surfaces as the adversarial example often being perceptually of a different class, although it originated from perturbing an image of a different class. Figure \ref{fig:op-abuse} demonstrates this failure of the $\ell_2$ norm, although there exists many more examples.

\begin{figure}[H]
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.5\linewidth]{images_op/digit_9.png}
  \caption{Unperturbed MNIST Data}
\end{subfigure}%
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.5\linewidth]{images_op/digit_9_adv.png}
  \caption{Overpowered Attack Data}
\end{subfigure}
\caption{Original and perturbed MNIST data under Overpowered Attack. $\ell_2$ difference of $3.8232$.}
\label{fig:op-abuse}
\end{figure}

The implications of the phenomenon of latent space attacks transmuting an data point of one class to another raises some issues about how we measure the effectiveness of latent space attacks, and whether they can be used at all. In terms of creating defenses for latent space attacks, it is difficult to gauge accuracy since it is hard to know what images have changed human perceptible classes. Similarly, it is difficult to perform adversarial training (either strengthening a classifier or a defense model) because this attack may introduce incorrectly classified images as true. One possible solution may be to constrain the $\ell_2$ distance further, or perhaps to use a different metric altogether. I leave this problem of latent space class transmutation as an open problem to the reader.

\section{Conclusions and Future Works} \label{conclusion}
In this paper, we propose the MAD-VAE model which improves upon the robustness of Defense-VAE through three new training loss functions. While applied here to Defense-VAE, our loss function and methodology is not limited to any model and can be applied to any generative model based defense algorithms. We empirically show that our methods are more effective against various attack mechanisms than Defense-VAE.

The only attack method which MAD-VAE does not improve upon the accuracy of Defense-VAE is the Overpowered \cite{jalal2017robust} Attack introduced by Jalal et al. While empirically our method fails to outperform Defense-VAE, in this paper we have raised concerns as to the applicability of the Overpowered Attack and other latent space attacks in their current state. It is possible that our drop in accuracy is due to an increase of images in which they have been transmuted to a different class.

An aspect that worth delving into would be adapting our training methods to other models such as Defense-GAN. Since it is training a GAN is not an easy task, it was not feasible to complete it within the limited time-frame given to do this research. Additionally, the choice of different hyperparameters is also worth tuning for a longer period of time. In order to achieve better performance, a larger training data set could be generated and could be used to better tune the hyperparameters. Additionally, additional research on latent space attacks could shed more light on avoiding the class transmutation problem.

\section*{Acknowledgement}
We are thankful for the help and insightful suggestions offered by Professor Shuyang Ling at New York University Shanghai.

\newpage

\bibliographystyle{unsrt}


\section{Introduction}
\label{sec:intro}
\vspace{-0.2cm}
Various interference such as environmental noise and room reverberation may result in heavy performance degradation in automatic speech recognition (ASR) and hearing assistance devices. Despite endless efforts in the last five decades, speech enhancement (SE), a type of technique to extract the clean version from the noisy speech, still remains a tough challenge under low signal-to-noise ratio (SNR) and nonstationary noise conditions~{\cite{loizou2013speech}}.

In recent years, the rapid development of deep neural networks (DNNs) has facilitated the research toward supervised speech enhancement algorithms~{\cite{wang2018supervised, kolbaek2020loss, xu2014regression}}. In the conventional supervised SE paradigm, DNN is leveraged to extract a clean feature from a noisy observation in the time-frequency (T-F) domain, $\emph{e.g.}$, ideal ratio mask (IRM)~{\cite{hummersone2014ideal}} or log-power spectrum (LPS)~{\cite{xu2014regression}}. With this aim, these methods usually focus on the estimation of the spectrum magnitude, and the noisy phase is left unaltered to reconstruct the waveform in the time domain.

More recently, the importance of phase begins to be emphasized as it is helpful to improve the perceptual quality under low SNR conditions~{\cite{paliwal2011importance}}. However, phase is usually difficult to estimate due to its unstructured characteristics. To this end, some SE algorithms are thus proposed and can be roughly divided into two categories. The first one tries to recover the phase information implicitly in the frequency domain, where both real and imaginary (RI) components of the spectrum are taken into optimization when using the complex spectral mapping based networks~{\cite{tan2019learning}}. The second one is to divert around direct phase estimation problem in the time domain, where the time-domain waveform is utilized as both the input and the output~{\cite{kolbaek2020loss, pandey2020densely}}. It has shown that noise and speech components tend to be more distinguishable in the T-F domain, and meanwhile, it is also easier for network to train in the T-F domain~{\cite{yin2020phasen}}, and thus this study focuses on SE in the T-F domain.

Multi-stage learning has been demonstrated to be more effective than single-stage methods in many tasks, like image deraining~{\cite{li2018recurrent}} and SE~{\cite{li2020speech, li2020recursive, hao2020masking, 9043689}}. In multi-stage learning protocol, the original difficult task is decomposed into multiple easier sub-tasks and the estimated target is progressively improved. Motivated by that, we propose a \textbf{C}omplex spectral mapping based \textbf{T}wo-\textbf{S}tage \textbf{N}etwork called CTS-Net for monaural speech enhancement. It is comprised of two sub-networks, namely \textit{coarse magnitude estimation network} (dubbed CME-Net) and \textit{complex spectrum refine network} (dubbed CSR-Net). In the first stage, the target spectral magnitude (TMS) is coarsely estimated by CME-Net, which is then coupled with noisy phase to obtain a coarsely estimated complex spectrum (dubbed coarse spectrum). In the second stage, the complex spectrum is further refined by CSR-Net, which recovers both the real and imaginary (RI) components. Note that, instead of explicitly estimating the complex spectrum, we only estimate the residual details in the second stage. This is because most noise interference has been removed in CME-Net and the second network aims to recover the clean speech phase, further suppress the residual noise, and restore some missing spectral details in the first stage.

The rationale behind such network design logic can be explained from two aspects. On one hand, it has been illustrated that an optimal solution for multi-task $\emph{w.r.t.}$ magnitude and phase can not be obtained simultaneously~{\cite{wang2020complex}}, especially in extremely low SNR condtions. For example, while the phase is continuously optimized by estimating RI components, magnitude estimation may gradually deviate its optimal optimization path. Therefore, the two-stage network topology can decouple the multi-task optimization problem w.r.t. magnitude and phase into two sub-tasks. On the other hand, it has been revealed that there exists a latent connection between magnitude and phase in the T-F domain~{\cite{griffin1984signal}}, and thus this relatively ``clean" magnitude obtained by the first stage can facilitate the phase recovery in the second stage.

The remainder of the paper is organized as follows. Section~{\ref{sec:problem-formulation}} formulates the problem. In Section~{\ref{sec:proposed-architecture}}, the proposed architecture is introduced in detail. Section~{\ref{sec:experiment-result}} gives the experimental results and analysis. Some conclusions are drawn in Section~{\ref{sec:conclusion}}.

\vspace{-0.4cm}
\section{PROBLEM FORMULATION}
\label{sec:problem-formulation}
\vspace{-0.2cm}
In the time domain, a single-microphone mixture $x$ is usually formulated as $x\left( t \right) = s\left( t \right) + n \left( t \right)$, where $s$ and $n$ denote clean and noise signal in the time index $t$. Taking the short-time Fourier transform (STFT) on both sides, we have:
\begin{equation}
\label{eqn:equa1}
X_{m, l} = S_{m, l} + N_{m, l},
\end{equation}
where $X_{m, l} = \lvert X_{m, l} \rvert e^{j\theta_{X_{m, l}}}$, $S_{m, l} = \lvert S_{m, l} \rvert e^{j\theta_{S_{m, l}}}$ and $N = \lvert N_{m, l} \rvert e^{j\theta_{N_{m, l}}}$ refer to the STFT representations of noisy, clean and noise, resepctively. $m$ and $l$ denote the frequency index and the time/frame index, respectively. For notation simplicity, we drop $\left(m, l\right)$ when no confusion arises.

Phase is usually left unchanged due to its unstructured characteristic for conventional SE algorithms. However, in Cartesian coordinates, phase can be implicitly represented by RI components, $\emph{e.g.}$, $\theta_{S} = \arctan \left(S_{i}/ S_{r}\right)$, which provides a promising direction for phase estimation. Recently, Tan $\emph{et.al.}$~{\cite{tan2019learning}} proposed a complex spectral mapping method with convolutional recurrent network (CRN), which took noisy and clean RI as inputs and targets, respectively. Assuming the network mapping function and its parameter set are $\mathcal{G}_{C}$ and $\phi_{C}$, respectively, the mapping process can be given by:
\begin{gather}
\label{eqn:equa2}
\left( \tilde{S}_{r}, \tilde{S}_{i} \right) = \mathcal{G}_{C}\left(X_{r}, X_{i};\phi_{C} \right),
\end{gather}
where $\left( \tilde{S}_{r}, \tilde{S}_{i} \right)$ denote the estimated RI.

\vspace{-0.2cm}
\section{PROPOSED ARCHITECTURE}
\label{sec:proposed-architecture}

\begin{figure}[t]
	\centering
	\centerline{\includegraphics[width=0.85\columnwidth]{mg_tcm.pdf}}
	\caption{Comparison between original TCM and proposed MG-TCM. (a) The diagram of original TCM. (b) The diagram of proposed MG-TCM.}
	\label{fig:gate-tcm}
	\vspace{-0.4cm}
\end{figure}
\vspace{-0.2cm}
\subsection{Modified gated temporal convolution module}
\label{sec:modified-gated}
\vspace{-0.0cm}
In~{\cite{luo2019conv}}, a convolution-based module named temporal convolution module (TCM) was proposed to replace traditional long-short term memory units (LSTMs) for better temporal sequence modeling. The diagram of TCM is shown in Fig.~{\ref{fig:gate-tcm}} (a). It consists of three major parts, namely input 1$\times$1-conv, dilated depthwise convolution (DD-conv) with kernel size 3, and output 1$\times$1-conv. Assuming the input size is $\left(256, T \right)$, where 256 and $T$ denote the channel and time axis, respectively. The number of channels is first doubled to 512, followled by a DD-conv. Then another 1$\times$1-conv layer is utilized to switch the channel back to 256. Residual connection is adopted to mitigate the gradient vanishing problem.

Despite the excellent performance of TCM, its drawback in heavy number of parameters is still noticeable. Specifically, although depthwise convolution is adopted to increase the parameter efficiency, as the channel dimension is increased to 512, the number of parameters in both input and output 1$\times$1-convs are still relatively large. To further ease the parameter burden, we propose a light-weight TCM variant named modified gated TCM (MG-TCM), as shown in Fig.~{\ref{fig:gate-tcm}} (b). Compared with original TCM (dubbed O-TCM), it has several modifications. Firstly, the channel dimension is decreased to 64 after the input 1$\times$1-conv, which notably decreases the parameter redundancy. Secondly, considering the limited ability for depthwise convolution in feature transform, we replace it with regular dilated convolution (D-conv) of kernel size 5 to grasp more information between adjacent frames. Thirdly, apart from the main branch, we also add another branch called gate branch, which has a similar structure as the main branch, except the sigmoid function is applied as the non-linearity to modulate the feature distribution. As shown in Fig.~{\ref{fig:gate-tcm}}, if neglecting the parameters for norm and activation layers, the number of parameters for O-TCM is 263,680, while only 53,248 for MG-TCM.

\begin{figure}[t]
	\centering
	\centerline{\includegraphics[width=0.85\columnwidth]{dmg_tcm.pdf}}
	\caption{Diagram of DMG-TCM. Norm and activation layers are omitted for convenience. ``CAT" denotes concatenation operation.}
	\label{fig:dmg-tcm}
	\vspace{-0.6cm}
\end{figure}

TCMs are usually stacked to grasp a larger temporal contextual field, where the dilation rate in each block is exponentially increased. If the dilation rate gets larger, then more attention is placed to long-term dependency, thus neglecting the correlation in the local regions. In this case, we further propose a dual version of MG-TCM named DMG-TCM to address this issue, as shown in Fig.~{\ref{fig:dmg-tcm}}. Compared with MG-TCM, it has two domains, namely primal and dual domain. The module detail within each domain is similar to that of MG-TCM except for the dilation rate. Suppose the dilation rate in the primal domain is $d_{\mathcal{P}} = 2^{r}$, then the dilation rate in the dual domain becomes $d_{\mathcal{D}} = 2^{M-r}$. In our method, $M = 5$ suffices empirically. When $d_{\mathcal{P}}$ in primal domain is small, $\emph{e.g.}$, $d_{\mathcal{P}} = 1$, the sequence in primal domain will pay more attention to local correlation, then $d_{\mathcal{D}}$ in dual domain becomes $32$, so more information among long-temporal correlation can be grasped. In this way, two domains complement with each other in the sequence modeling process.

\vspace{-0.4cm}
\subsection{Smoothed dilated convolution}
\label{sec:sd-conv}
Despite the effectiveness of dilated convolution, it may have the so-called ``gridding artifacts"~{\cite{hamaguchi2018effective}}. That is, consecutive units in the output are computed from completely separate sets of units in the input, and the actual receptive field is different. To handle this problem,Wang and Ji $\emph{et.al.}$~{\cite{wang2018smoothed}} proposed the smoothed dilated convolutions (SD-convs) where the interaction was added among adjacent input units before adopting dilation convolution. A separable and shared (SS) convolution is applied before calculating dilated convolution. Here ``separable" means that each channel deals with its kernel filters independently, and ``shared" denotes that the parameters of kernel filters are shared for different input and output channels. In this paper, all the D-convs are replaced by SD-convs to alleviate the artifacts.

\begin{figure}[t]
	\centering
	\centerline{\includegraphics[width=0.92\columnwidth]{two_stage.pdf}}
	\caption{The diagram of proposed two-stage algorithm. (a) The processing flow of the algorithm. (b) The model topology of CME-Net. (c) The model topology of CSR-Net.}
	\label{fig:two-satge}
	\vspace{-0.6cm}
\end{figure}
\vspace{-0.4cm}
\subsection{Two-stage approach}
\label{sec:two-stage}
Our proposed two-stage algorithm is illustrated in Fig.~{\ref{fig:two-satge}}. It includes two processing stages. In the first stage, CME-Net receives the magnitude of noisy spectrum to coarsely estimate the magnitude of the clean speech, which is then coupled with noisy phase to obtain a coarse complex spectrum. In the second stage, the real and imaginary components for both coarse and noisy spectra are concatenated as the inputs of CSR-Net. Instead of directly refining the complex spectrum, it only captures the spectral details which might be lost in the first stage. Formally, the calculation process is formulated as:
\vspace{-0.2cm}
\begin{gather}
\label{eqn:equa3}
\lvert \tilde{S}^{cm} \rvert = \mathcal{G}_{cm} \left( \lvert X \rvert; \phi_{cm} \right),\\
\tilde{S}^{cm}_{r} = \Re \left(\lvert S^{cm} \rvert e^{j\theta_{X}} \right),
\tilde{S}^{cm}_{i} = \Im \left(\lvert S^{cm} \rvert e^{j\theta_{X}} \right),\\
\left( \tilde{S}^{cs}_{r}, \tilde{S}^{cs}_{i} \right) = \left( \tilde{S}^{cm}_{r}, \tilde{S}^{cm}_{i} \right) + \mathcal{G}_{cs} \left(\tilde{S}^{cm}_{r}, \tilde{S}^{cm}_{i}, X_{r}, X_{i}; \phi_{cs} \right),
\end{gather}
where $\mathcal{G}_{cm}$ and $\mathcal{G}_{cs}$ denote the mapping functions for CME-Net and CSR-Net, respectively. $\phi_{cm}$ and $\phi_{cs}$ denote the parameter sets for the first and second network, respectively. $\tilde S^{cm}$ and $\tilde S^{cs}$ are the estimated outputs from CME-Net and CSR-Net, respectively. $\Re$ and $\Im$ refer to real and imaginary operations, respectively.

The network details are shown in Fig.~{\ref{fig:two-satge}}(b)-(c). The overall topology is similar to that of~{\cite{tan2019learning}}, which includes three major components, namely convolutional encoder, decoder, and sequence modeling module. Instead of using LSTMs as the sequence module, we adopt TCMs for better sequence learning. For CME-Net, we stack 18 MG-TCMs while 12 DMG-TCMS are adopted in CSR-Net~{\footnote{In this paper, 6 TCMs form a large group, where the dialtion rate in each group is (1, 2, 4, 8, 16, 32)}. This is because DMG-TCM is more powerful to capture both long and short-term temporal information than MG-TCM. In both encoders and decodes, 5 convolutional blocks are adopted, each of which includes one (de)convolution layer, instance normalization~{\cite{ulyanov2016instance}}, and Parametric ReLU (PReLU)~{\cite{he2015delving}}. The kernel size for each (de)convolution layers is (2, 3) except (2, 5) for the first layer in the time and frequency axis, respectively. The stride is set to (1, 2) in the time and frequency axis. The number of channels for each intermediate layer is 64. In the first network, Softplus~{\cite{glorot2011deep}} is used as the output activation function to obtain the magnitude of the spectrum and linear function is used to obtain the RI in the second stage. Note that similar to~{\cite{tan2019learning}}, two decoders are used in CSR-Net to obtain RI estimation.
\vspace{-0.2cm}
\subsection{Loss function}
\label{loss-function}
We take the following strategy to train the network. Firstly, we separately train CME-Net until convergence, and the loss is defined as:
\begin{gather}
\label{eqn:equa4}
\mathcal{L}_{cm} =\left \| \lvert \tilde S^{cm} \rvert - \lvert S \rvert \right \|^{2}_{F},
\end{gather}
Afterward, both the first and the second stages are jointly trained, where the first stage is initialized with the pretrained model, and the overall loss is defined as:
\begin{gather}
\label{eqn:equa5}
\mathcal{L} = \mathcal{L}_{cs}^{RI} + \mathcal{L}_{cs}^{Mag} + \lambda \mathcal{L}_{cm},\\
\mathcal{L}_{cs}^{RI} = \left \| \tilde{S}^{cs}_{r} - S_{r} \right \|^{2}_{F} + \left \| \tilde{S}^{cs}_{i} - S_{i} \right \|^{2}_{F},\\
\mathcal{L}_{cs}^{Mag} = \left \| \sqrt{ \lvert \tilde{S}_{r}^{cs}\rvert ^2 + \lvert \tilde{S}_{i}^{cs} \rvert ^2  } - \sqrt{ \lvert S_{r}\rvert ^2 + \lvert S_{i} \rvert ^2} \right \|_{F}^{2},
\end{gather}
where $\mathcal{L}^{RI}$ and $\mathcal{L}^{Mag}$ denote the loss optimization toward RI and that toward magnitude, respectively. $\lambda \in [0,1]$ controls the loss weight of the first network. In this paper, $\lambda$ is set to 0.1.

\vspace{-0.2cm}
\section{EXPERIMENTAL RESULTS}
\label{sec:experiment-result}
\vspace{-0.2cm}

\subsection{Datasets}
We use the WSJ0-SI84 dataset~{\cite{paul1992design}} for evaluation, which includes 7138 utterances by 83 speakers (42 males and 41 females). 5428, and 957 utterances with 77 speakers are split for training and validation, respectively. Two types are set for test. For the first type, the speaker information is within the training dataset (dubbed seen speaker), while for the second type, the speaker information is untrained (dubbed unseen speaker). Each of the types includes 150 utterances. We randomly select 20,000 noises from the DNS-Challenge~{\footnote{https://github.com/microsoft/DNS-Challenge}} to obtain a 55 hours noise set for training. During each mixed process, a random cut is generated to obtain a noise vector, which is subsequently mixed with randomly selected clean utterance. The SNR range for training is [-5\rm{dB}, 0\rm{dB}] with the interval 1\rm{dB}. As a result, totally 50,000, 4000 noisy-clean pairs are established for training and validation, respectively. The total duration for training set is about 100 hours.

For test, we select two noises from NOISEX92~{\cite{varga1993assessment}}, namely babble and factory1. Three SNRs are used for model evaluation, namely -5\rm{dB}, 0\rm{dB} and 5\rm{dB}.

\renewcommand\arraystretch{1.62}
\begin{table*}[t]
	\setcounter{table}{1}
	\caption{Objective result comparisons among different models in terms of PESQ, ESTOI and SDR for both seen and unseen speaker cases. \textbf{BOLD} indicats the best score in each case.}
	\Huge
	\centering
	\resizebox{0.82\textwidth}{!}{
		\begin{tabular}{cc|c|cccccccc|cccccccc|cccccccc}
			\hline
			 &\textbf{Metrics} &\multirow{3}*{\rotatebox{90}{\textbf{Causality}}} &\multicolumn{8}{c|}{\textbf{PESQ}} &\multicolumn{8}{c|}{\textbf{ESTOI(\%)}} &\multicolumn{8}{c}{\textbf{SDR(dB)}} \\
			\cline{1-2}\cline{4-27}
			& \textbf{Noises} &  &\multicolumn{4}{c|}{\textbf{Babble}} &\multicolumn{4}{c|}{\textbf{Factory1}} &\multicolumn{4}{c|}{\textbf{Babble}} &\multicolumn{4}{c|}{\textbf{Factory1}} &\multicolumn{4}{c|}{\textbf{Babble}} &\multicolumn{4}{c}{\textbf{Factory1}} \\
			\cline{1-2}\cline{4-27}
			& \textbf{SNR(dB)} & &-5 &0 &5 &\multicolumn{1}{c|}{Avg.} &-5 &0 &5 &Avg. &-5 &0 &5 &\multicolumn{1}{c|}{Avg.} &-5 &0 &5 &\multicolumn{1}{c|}{Avg.} &-5 &0 &5 &\multicolumn{1}{c|}{Avg.} &-5 &0 &5 &Avg.\\
			\cline{1-27}
			\multirow{6}*{\rotatebox{90}{\textbf{Seen speaker}}}
			&\multicolumn{1}{|c|}{\textbf{Noisy}} &-  &1.64 &1.88 &2.19 &\multicolumn{1}{c|}{1.90} &1.51 &1.79 &2.13 &1.81 &29.32 &42.96 &58.74 &\multicolumn{1}{c|}{43.67} &28.89 &44.41 &60.73 &44.68 &-4.91 &0.04 &5.03 &\multicolumn{1}{c|}{0.05} &-4.90 &0.04 &5.03 &0.06\\
			\cline{2-27}
			&\multicolumn{1}{|c|}{\textbf{CRN}} &$\checkmark$ &1.86 &2.34 &2.76 &\multicolumn{1}{c|}{2.32} &2.01 &2.45 &2.84 &2.43 &49.10 &65.38 &77.35 &\multicolumn{1}{c|}{63.94} &48.63 &65.99 &77.87 &64.16 &1.55 &6.28 &10.45 &\multicolumn{1}{c|}{6.09} &3.36 &7.51 &11.18 &7.35\\
			\cline{2-27}
			&\multicolumn{1}{|c|}{\textbf{TCNN}} &$\checkmark$  &1.91 &2.42 &2.83 &\multicolumn{1}{c|}{2.39} &2.06 &2.50 &2.84 &2.47 &55.12 &70.91 &81.13 &\multicolumn{1}{c|}{69.05} &54.70 &71.04 &80.65 &68.80 &4.55 &9.33 &12.98 &\multicolumn{1}{c|}{8.95} &6.07 &10.07 &13.15 &9.76\\
			\cline{2-27}
			&\multicolumn{1}{|c|}{\textbf{GCRN}} &$\checkmark$  &2.02 &2.53 &2.91 &\multicolumn{1}{c|}{2.49} &2.14 &2.63 &3.00 &2.59 &54.18 &70.07 &80.37 &\multicolumn{1}{c|}{68.21} &54.14 &71.41 &81.12 &68.89 &3.96 &8.44 &11.92 &\multicolumn{1}{c|}{8.11} &6.00 &9.57 &12.66 &9.41\\
			\cline{2-27}
			&\multicolumn{1}{|c|}{\textbf{CME-Net(Pro.)}} &$\checkmark$ &1.98 &2.47 &2.86 &\multicolumn{1}{c|}{2.44} &2.12 &2.55 &2.91 &2.53 &51.77 &67.46 &78.72 &\multicolumn{1}{c|}{65.99} &51.44 &67.69 &78.76 &65.96 &2.32 &6.85 &10.87 &\multicolumn{1}{c|}{6.68} &3.76 &7.76 &11.47 &7.66\\
			\cline{2-27}
			&\multicolumn{1}{|c|}{\textbf{CTS-Net(Pro.)}} &$\checkmark$  &\textbf{2.18} &\textbf{2.75} &\textbf{3.13} &\multicolumn{1}{c|}{\textbf{2.69}} &\textbf{2.31} &\textbf{2.81} &\textbf{3.14} &\textbf{2.76} &\textbf{60.10} &\textbf{75.37} &\textbf{84.12} &\multicolumn{1}{c|}{\textbf{73.19}} &\textbf{59.82} &\textbf{75.05} &\textbf{83.37} &\textbf{72.74} &\textbf{5.46} &\textbf{10.21} &\textbf{13.63} &\multicolumn{1}{c|}{\textbf{9.77}} &\textbf{6.93} &\textbf{10.74} &\textbf{13.70} &\textbf{10.46}\\
			\hline
			\multirow{6}*{\rotatebox{90}{\textbf{Unseen speaker}}}
			&\multicolumn{1}{|c|}{\textbf{Noisy}} &$\checkmark$ &1.55 &1.83 &2.15 &\multicolumn{1}{c|}{1.84} &1.43 &1.74 &2.10 &1.76 &26.60 &39.60 &54.43 &\multicolumn{1}{c|}{40.21} &26.66 &41.01 &57.03 &41.57 &-4.91 &-0.05 &5.03 &\multicolumn{1}{c|}{0.02} &-4.93 &0.04 &5.03 &0.05\\
			\cline{2-27}
			&\multicolumn{1}{|c|}{\textbf{CRN}} &$\checkmark$ &1.76 &2.28 &2.70 &\multicolumn{1}{c|}{2.24} &1.95 &2.40 &2.78 &2.38 &45.31 &62.18 &75.12 &\multicolumn{1}{c|}{60.87} &46.03 &63.07 &75.95 &61.68 &1.43 &6.37 &10.63 &\multicolumn{1}{c|}{6.14} &3.37 &7.63 &11.45 &7.48\\
			\cline{2-27}
			&\multicolumn{1}{|c|}{\textbf{TCNN}} &$\checkmark$  &1.83 &2.36 &2.78 &\multicolumn{1}{c|}{2.32} &2.00 &2.45 &2.80 &2.41 &51.88 &68.92 &79.86 &\multicolumn{1}{c|}{66.89} &51.59 &68.95 &79.55 &66.70 &4.33 &9.35 &13.10 &\multicolumn{1}{c|}{8.93} &6.16 &10.13 &13.39 &9.89\\
			\cline{2-27}
			&\multicolumn{1}{|c|}{\textbf{GCRN}} &$\checkmark$ &1.94 &2.49 &2.90 &\multicolumn{1}{c|}{2.44} &2.13 &2.61 &2.98 &2.57 &52.19 &69.09 &79.88 &\multicolumn{1}{c|}{67.05} &52.85 &70.12 &80.30 &67.76 &3.97 &8.73 &12.37 &\multicolumn{1}{c|}{8.36} &5.85 &9.81 &13.03 &9.56\\
			\cline{2-27}
			&\multicolumn{1}{|c|}{\textbf{CME-Net(Pro.)}} &$\checkmark$  &1.91 &2.42 &2.81 &\multicolumn{1}{c|}{2.38} &2.08 &2.49 &2.86 &2.48 &48.53 &65.20 &76.84 &\multicolumn{1}{c|}{63.53} &49.07 &65.22 &77.23 &63.84 &2.47 &7.22 &11.19 &\multicolumn{1}{c|}{6.96} &4.04 &7.97 &11.81 &7.94\\
			\cline{2-27}
			&\multicolumn{1}{|c|}{\textbf{CTS-Net(Pro.)}} &$\checkmark$ &\textbf{2.11} &\textbf{2.72} &\textbf{3.11} &\multicolumn{1}{c|}{\textbf{2.65}} &\textbf{2.29} &\textbf{2.77} &\textbf{3.12} &\textbf{2.72} &\textbf{57.82} &\textbf{74.14} &\textbf{83.15} &\multicolumn{1}{c|}{\textbf{71.70}} &\textbf{58.07} &\textbf{73.59} &\textbf{82.70} &\textbf{71.45} &\textbf{5.40} &\textbf{10.46} &\textbf{13.96} &\multicolumn{1}{c|}{\textbf{9.94}} &\textbf{7.24} &\textbf{10.96} &\textbf{14.14} &\textbf{10.78}\\
			\hline
	\end{tabular}}
	\label{tbl:objective-results}
	\vspace{-0.4cm}
\end{table*}
\vspace{-0.2cm}
\subsection{Parameter setup}
\label{parameter-setup}
All the utterances are sampled at 16kHz. The 20ms Hanning window is utilized, with 50\% overlap between adjacent frames. 320 point FFT is used. Both models are optimized by Adam~{\cite{kingma2014adam}} with $\beta_{1} = 0.9, \beta_{2} = 0.999$. In the first stage, the learning rate (LR) is set to 0.001. In the second stage, the pre-trained model in the first stage is fine-tuned with LR = 0.0001, while 0.001 for the second model. The batch size is set to 16 at an utterance level, where the maximum utterance length is chunked to 8 seconds for training stability.
\vspace{-0.2cm}
\subsection{Baselines}
We adopt three state-of-the-art baselines for comparison, namely CRN~{\cite{tan2018convolutional}}, TCNN~{\cite{pandey2019tcnn}} and GCRN~{\cite{tan2019learning}}. CRN is a typical convolutional recurrent network with encoder-decoder architecture, and only magnitude is estimated. GCRN is an advanced complex spectral mapping network based on CRN, where both magnitude and phase are estimated. Note that both RI and magnitude are optimized with $\mathcal{L}^{RI}$ and $\mathcal{L}^{Mag}$ for GCRN. For TCNN, the waveform is directly used as both input and target, and stacked TCMs are adopted for sequence modeling. All the models are trained with causal configurations, $\emph{i.e.}$, no future information is involved. The sound demo samples are available online~{\footnote{https://github.com/Andong-Li-speech/CTS-Net}}.

\vspace{-0.2cm}

\subsection{Results and analysis}
\label{results-analysis}
\vspace{-0.2cm}
In this study, PESQ~{\cite{rix2001perceptual}}, ESTOI~{\cite{jensen2016algorithm}} and SDR~{\cite{vincent2007first}} are chosen as objective measurement metrics.
\vspace{-0.0cm}

\renewcommand\arraystretch{0.6}
\begin{table}[t]
	\setcounter{table}{0}
	\caption{Ablation study $\emph{w.r.t.}$ TCM type and SD-conv. TCM(I) and TCM(II) denote the TCMs used in the first and second stage, respectively. ``O", ``MG" and ``DMG" denote O-TCM, MG-TCM and DMG-TCM. All the values are averaged on the test set.}
	\large
	\centering
	\resizebox{0.40\textwidth}{!}{
		\begin{tabular}{ccccccc}
			\toprule
			Models
			& TCM(I)  &TCM(II) &SD  &PESQ &ESTOI(\%) &SDR (dB)\\
			\midrule
			Noisy &- &- &- &1.83 &42.53 &0.30\\
			\midrule
			CME-Net  &O &- &$\times$  &2.40 &64.03 &7.38 \\
			\midrule
			CME-Net &MG &- &$\times$ &2.44 &64.53 &7.35\\
			\midrule
			CTS-Net &MG &MG &$\checkmark$ &2.70 &71.89 &10.16 \\
			\midrule
			CTS-Net &MG &DMG &$\checkmark$ &$\textbf{2.71}$ &$\textbf{72.27}$ &$\textbf{10.24}$\\
			\midrule
			CTS-Net &MG &DMG &$\times$ &2.68 &71.75 &10.11\\
			\bottomrule
		\end{tabular}}
	\label{tbl:ablation-study}
\vspace{-0.4cm}
\end{table}
\vspace{-0.4cm}

\subsubsection{Ablation study}
\label{ablation-study}
\vspace{-0.0cm}
We investigate the effects of different TCMs and SD-conv, as shown in Table~{\ref{tbl:ablation-study}}. From the results, we can have the following observations. Firstly, MG-TCM achieves overall better performance than O-TCM. For example, when only the first stage is trained, CME+MG-TCM obtains 0.04 and 0.50\% improvements in PESQ and ESTOI than CME+O-TCM, which validates the parameter redundancy for original TCM. Secondly, when MG-TCM is replaced by its dual version in the second stage, consistently better performance is achieved in all three metrics. Thirdly, compared with regular dilated convolution, SD-conv provides 0.03, 0.52\% and 0.13\rm{dB} improvements in PESQ, ESTOI and SDR, respectively.
\vspace{-0.0cm}
\subsubsection{Comparison with baselines}
\label{comparison-with-baseline}
\vspace{-0.2cm}
The results of different models are shown in Table~{\ref{tbl:objective-results}}. Note that CME-Net and CTS-Net are the best configurations from the ablation study. One can have the following observations. Firstly, CME-Net significantly outperforms CRN in different cases. For example, for seen speaker, CME-Net provides average 0.12 and 0.10 PESQ improvements than CRN on babble and factory1 noises, while 2.05\% and 1.80\% improvements in ESTOI. This indicates the superior performance of MG-TCMs at sequence modeling capability than naive LSTMs. Secondly, compared with CME-Net, when the second stage is leveraged to refine the spectrum, large metric improvements are achieved. For example, going from CME-Net to CTS-Net, average 0.25 and 6.94\% improvements are achieved in PESQ and ESTOI, respectively. This reveals the necessity and significance of the second stage in improving the speech quality and intelligibility. Thirdly, the proposed two-stage model consistently surpasses all the baselines. For example, compared with GCRN, average 0.18, 4.29\%, and 1.38\rm{dB} metric improvements are obtained in terms of PESQ, ESTOI and SDR, respectively.

Table~{\ref{tbl:model-parameters}} summarizes the number of trainable parameters among different models. One can find that the CME-Net has the smallest number of parameters among different models. As CTS-Net includes two cascaded sub-networks, the number of parameters is relatively large. Nevertheless, it is still less than another three baselines.

\renewcommand\arraystretch{1.0}
\begin{table}[t]
	\setcounter{table}{2}
	\caption{The number of trainable parameters among different models. The unit is million. $\textbf{BOLD}$ denotes the lowest trainable parameters. }
	\centering
	\footnotesize
	\resizebox{0.45\textwidth}{!}{
	\begin{tabular}{c|c|c|c|c|c}
		\hline
		Model &CRN &GCRN &TCNN &CME-Net &CTS-Net\\
		\hline
		Para. (million) &17.59 &9.06 &5.08 &\textbf{1.96} &4.99 \\
		\hline
	\end{tabular}}
	\label{tbl:model-parameters}
	\vspace*{-\baselineskip}
	\vspace{-0.1cm}
\end{table}

\vspace{-0.2cm}
\section{CONCLUSIONS}
\label{sec:conclusion}
\vspace{-0.2cm}
In this work, we propose a two-stage algorithm for monaural noise reduction in the complex domain. In the first stage, the magnitude of spectrum is estimated, which is coupled with noisy phase to obtain a coarse complex spectrum. In the second stage, the spectral details are captured, which further refines the magnitude and phase information simultaneously. Additionally, a modified TCM is proposed, which can achieve better performance than previous counterparts in sequence learning while decreasing the parameter redundancy by a large margin. Experimental results indicate that the proposed algorithm consistently outperforms previous powerful methods while still enjoying a small parameter capacity.
\vspace{-0.2cm}
\section{ACKNOWLEDGMENT}
\label{sec:acknowledge}
\vspace{-0.2cm}
This work was supported by National Key R\&D Program of China. This work was also supported by IACAS Young Elite Researcher Project under no.QNYC201813. We would like to thank Cunhang Fan at Institute of Automation, and Shan You at SenseTime Research for constructive comments.

\vfill\pagebreak

\bibliographystyle{IEEEbib}


\section*{Method}
\noindent\textbf{Galaxy cluster sample selection.}
We construct our sample of high-$z$ galaxy clusters using the latest {\it Planck} Sunyaev-Zel'dovich (SZ) catalog, i.e. PSZ2 \cite{planckcoll16}. This catalog provides a reliable cluster mass estimation for galaxy clusters up to $z\sim1$. Although the exact effect of selection biases is still being debated and may overestimate the observed fraction of clusters with radio halos,
there is no clear consensus in the literature that suggests the SZ sample chosen here disproportionately favours merging over relaxed clusters \cite{eckert+11,rossetti+17,andrade-santos+2017}. We select all the objects at $\rm DEC\geq 20^\circ$, to match the part of the sky with the best LOFAR sensitivity, and $z\geq0.6$. No mass-threshold has been applied in our cluster selection.
These selection criteria result in a sample of 30 objects (see Supplementary Table 1). The sample includes clusters where precise redshift measurements became available recently \cite{amodeo+18,barrena+18,burenin+18,sereno+18,streblyanska+18,vanderburg+16,zohren+19} and a previous-discovered {\it Planck} cluster, i.e. PLCK\,G147.3--16.6 \cite{vanweeren+14}. We used the optical images from the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS \cite{panstarss16}) to refine the accuracy of the cluster center coordinates, using the brightest cluster galaxy (BCG), when detected, or approximately the center of the distribution of the galaxies. We also inspected the {\it Chandra} and/or {\it XMM-Newton} X-ray image, when available, to verify the correspondence of the cluster coordinates with the center of the thermal gas distribution. Among these 30 {\it Planck} clusters, 21 objects were already covered as part of the LOFAR Two-metre Sky Survey (LoTSS \cite{shimwell+19}). From these 21 clusters one was affected by bad ionospheric observing conditions (i.e. PSZ2\,G088.98+55.07) and one lays on the same line of sight of a cluster at $z=0.3$ (i.e. PSZ2\,G097.52+51.70). Those are therefore excluded from the final sample (see Extended Data Figure 1), consisting of 19 objects at median redshift of $z=0.7$.

\bigskip
\bigskip

\noindent\textbf{LOFAR observations, data reduction, images and flux measurements.}
The LOFAR observations were carried out together with the LOFAR Two-metre Sky Survey (LoTSS \cite{shimwell+19}) in the 120--168 MHz frequency range. The survey consists in 8 hours of observation for each pointing, with a field of view of $\sim2.6^\circ$, full resolution of about $6''$, and median noise of about $70~\mu$Jy beam$^{-1}$. Given the large sky coverage already achieved by the survey, our targets were often observed by more than one pointing (see Supplementary Table 2), which further improved our signal to noise. For each LoTSS pointing, we performed standard data reduction \cite{shimwell+19}, which includes direction-independent and dependent calibration, and imaging of the full LOFAR field of view using \texttt{prefactor} \cite{vanweeren+16,williams+16,degasperin+19}, \texttt{killMS} \cite{tasse14,smirnov+tasse15} and \texttt{DDFacet} \cite{tasse+18}. We additionally improved the quality of the calibration, using the products of the pipeline, subtracting all the sources outside a region of $\sim15^\prime\times15^\prime$ surrounding the target, and performing extra phase and amplitude self-calibration loops in each target sub-field. To verify the quality of the calibration, each pointing is imaged independently using \texttt{WSClean} \cite{offringa+14}, with the wideband deconvolution mode (\texttt{channelsout=6}). For all our imaging, we employ automatic clean masks for the deconvolution. We utilize a $3\sigma_{\rm rms}$ masking threshold, with $\sigma_{\rm rms}$ the local map noise, and cleaning down to the $1\sigma_{\rm rms}$ level inside the mask. The images have a central frequency of 144 MHz, and are provided in the Extended Data Figure 1.

\bigskip

\noindent All the clusters in our sample were carefully visually inspected in the full-resolution image to search for extended radio emission.
To emphasize the presence of diffuse radio emission, we also produced low-resolution source-subtracted images. We first applied an $uv$-cut to the data, to filter out emission associated to linear sizes larger than 500 kiloparsecat the cluster redshift and
to create a clean component model of the compact sources. For PSZ2\,G089.39+69.36 we apply an $uv$-cut corresponding to 400 kpc.
During this step, we employed multiscale deconvolution \cite{offringa+17}, using scales of $\rm [0,4,8,16]\times pixelscale$ (with the pixel size of $1.5''$) to properly include and subtract somewhat extended radio galaxies. In addition, for the automatic deconvolution we lower the mask threshold to $1\sigma_{\rm rms}$ to subtract the faintest contaminating sources (see Supplementary Figure 1).
Finally, we subtracted the compact sources model from the visibilities
and tapered the $uv$-data to $\sim15''$ resolution. It is important to note that in case of an extended radio galaxy, with linear size $\gtrsim 500$~kpc, this method cannot properly subtract the radio emission from the $uv$-data. For that reason, the extended double-lobed radio galaxy just north of PLCKG147.3--16.6 has been manually excluded and blanked in the low-resolution source-subtracted image (Figure 1 and Extended Data Figure 1).
\bigskip

\noindent We used these low-resolution source-subtracted images to define the largest linear size (LLS) of the diffuse radio emission, following the $2\sigma_{\rm rms}$ contour. Systems with $\rm LLS<500~kpc$ were not included in our statistics.
Roundish, centrally-located significantly detected sources (see the next paragraph) that follow the ICM distribution were classified as radio halos; elongated, peripheral structures were classified as radio relics. We did not identify clear examples of fossil radio plasma sources (i.e. ``radio phoenices''), these typically have small sizes (i.e. 200~kpc or smaller), irregular shapes, and are not centered on the cluster center \cite{degasperin+17,mandal+20}. Similarly, we did not find clear examples of  ``mini-halos'', as these diffuse sources usually have LLS smaller than 500 kpc and are located in dynamically relaxed systems \cite{giacintucci+17}. Uncertain classifications,
due to contamination from extended radio galaxies in the cluster region and/or  ambiguous shapes,
were excluded from our final statistical analysis.
If we measure the LLS following the $3\sigma_{\rm rms}$ contour, we would classify 8 radio halos instead of 9 (see Supplementary Table 3).

\bigskip

\noindent To determine the integrated flux densities for the radio halos, we first measured the total cluster radio flux densities (i.e. compact sources and diffuse emission) from the low-resolution image, following the $2\sigma_{\rm rms}$ radio contour to fully cover the extension of the diffuse radio emission (see Extended Data Figure 1). In the Supplementary Table 3, we also report the integrated flux density measurements following the $3\sigma_{\rm rms}$ radio contour. %
We then mathematically subtracted the flux densities of the compact sources within the cluster region, measured in the full-resolution $uv$-cut image, to obtain the flux density associated to the diffuse radio sources \cite{cassano+19}.
For PSZ2\,G091.83+21.16, which also hosts a bright radio relic, we visually separated it from the radio halo region (see Supplementary Figure 2).
The uncertainties on the flux densities are estimated adding in quadrature the $\sigma_{\rm rms}$ statistical uncertainty, the 15\% systematic error associated to the LOFAR flux scale calibration, $f_{\rm cal}$ \cite{shimwell+19}, and the uncertainty due to the source subtraction in the cluster region, $\sigma_{\rm sub}$ \cite{cassano+13}, according the following Equation:

\begin{equation}
\Delta S = \sqrt{(f_{\rm cal} S_{\rm 144MHz})^2 + \sigma_{\rm rms}^2 N_{\rm beam,h} + \sigma_{\rm sub}^2} \, .
\end{equation}
Here, $\sigma_{\rm sub}^2=\sum_i\sigma_{\rm rms}^2 N_{{\rm beam,s}_i}$, and $N_{\rm beam, h}$ and  $N_{{\rm beam, s}_i}$ are the number of beams covering the halo and the subtracted $i$ sources, respectively. Uncertainties associated to possibly missed faint compact sources are not included, but these should be smaller than the systematic error associated to the LOFAR flux scale calibration.

\bigskip

\noindent We confirmed the presence of diffuse radio emission when the measured integrated radio flux density is larger than five times the flux uncertainty, i.e. $\frac{S_{\rm 144MHz}}{\Delta S}\geq 5$ (see Column 5 in Table 1).
Flux density measurements directly obtained from the low-resolution source-subtracted images are also provided in the Supplementary.

\bigskip

\noindent Finally, we calculated the $k$-corrected radio luminosity at 1.4 GHz:

\begin{equation}
P_{\rm 1.4GHz} =  4\pi D_L^2 \dfrac{S_{\rm 1.4GHz}}{(1+z)^{\alpha+1}} \, ,
\end{equation}
with $D_L$ the luminosity distance assuming standard $\Lambda$CDM cosmology (i.e. $H_{0} = 70$ km s$^{-1}$~Mpc$^{-1}$, $\Omega_{m} = 0.3$, and $\Omega_{\Lambda} = 0.7$) and $S_{\rm 1.4GHz}=S_{\rm 144MHz}\left (\frac{\rm 1.4GHz}{\rm 144MHz} \right )^\alpha$. Since the spectral indices of the radio halos are not measured, we assume a Gaussian distribution of
spectral indices $\alpha$ with mean value $-1.5$ and scatter $0.3$ for computing the radio luminosity. These values were chosen to encompass the typical range of spectral indices found in radio halos \cite{vanweeren+19}, and to take into account that they could be steeper than local systems \cite{cassano+06}. Only for PLCK\,G147.3--16.6, we used the spectral index obtained combining our observations with the GMRT literature data \cite{vanweeren+14} (see Supplementary). The radio luminosities are listed in Table 1,
whose uncertainties are obtained by performing Monte Carlo simulations including the flux density uncertainties and the spectral index uncertainty.

\bigskip
\bigskip

\noindent\textbf{X-ray observations and data reduction.}
Twelve clusters among our LOFAR sample have also {\it Chandra} and/or {\it XMM-Newton} observations (see Supplementary Table 4 and Extended Data Figure 2).

\bigskip

\noindent {\it Chandra} data were reduced using the \texttt{chav} software package (\url{http://hea-www.harvard.edu/\textasciitilde alexey/CHAV/}) with \texttt{CIAO v4.6} and applying the \texttt{CALDB v4.7.6} calibration files \cite{vikhlinin+05}. The processing includes the application of gain maps to calibrate photon energies, filtering out counts with ASCA grade 1, 5, or 7 and bad pixels, and a correction for the position-dependent charge transfer inefficiency (CTI). Periods with count rates with a factor of 1.2 above and 0.8 below the mean count rate in the 6--12 keV band were also removed. Standard blank-sky files were used for background subtraction. The final exposure-corrected images were made in the 0.5--2.0~keV band using a pixel binning of a factor of 4 (i.e $2^{\prime\prime}$) and combining the different ObsIDs together where relevant (i.e., PSZ2\,G160.83+81.66).

\bigskip

\noindent {\it XMM-Newton} observations were reduced using version 17.0.0 of the dedicated Science Analysis System (SAS). After converting the observation data files to unfiltered event lists, we extracted light curves in bins of 100 seconds in the energy ranges 10--12~keV for the MOS and 10--14~keV for the pn detectors onboard. Good time intervals were created excluding periods when the mean count rates in these light curves were different from the corresponding mean by more than 2 standard deviations. Images were extracted in the 0.4--7.0 keV band, using only the highest quality (FLAG==0) single to quadruple MOS events (PATTERN$\leq$12) and single to double pn events (PATTERN$\leq$4). Images were weighted by the respective exposure maps for each detector and combined, accounting for a factor $\sim$2.5 difference between the typical expected count rates for MOS and pn at a given source flux (due to the slightly different effective areas and the fact that the light focused by two of the XMM-Newton telescopes is split evenly between the MOS and RGS detectors). Since we are interested mostly in the morphology of high-redshift, thus relatively compact, clusters positioned close to the aim point, the effects of vignetting and instrumental particle background are minimal and have not been corrected for this analysis.

\bigskip
\bigskip

\noindent\textbf{Turbulent energy flux at high redshift.}
According to re-acceleration models, relativistic electrons in radio halos are re-accelerated by turbulence. Assuming that turbulence in the ICM is mainly injected by cluster mergers, then the turbulence injection rate is $\rho v_i^3/L_{\rm inj}$ \cite{casano+brunetti05}, with $v_i$ the impact velocity between the two merging clusters, $\rho$ the cluster mean density and $L_{\rm inj}$ the injection scale of turbulence (which can be assumed to be of the same order of the cluster size, about 1 Mpc at both $z=0.2$ and $z=0.7$).
The relative impact velocity of two subclusters with virial masses $M_v$ and $\Delta M_v$ which collide (at a distance $R_{v}$ between the centers) starting from an initial distance $d_0$ with zero velocity is given by \cite{sarazin02}:

\begin{equation}
v_i\simeq \sqrt{2G\frac{(M_{v}+\Delta M_v)}{R_{v}}
\left (1-\frac{1}{\eta_{v}}\right)} \, .
\label{eq:vi}
\end{equation}
Here, $d_0=\eta_v R_{v}$, $\eta_v\simeq 4\big(\frac{ M_{v}+ \Delta M} {M_{v}}\big)^{1/3}$, and $R_{v}$ is the virial radius of the main cluster, i.e., the radius at which the ratio between the average density in the cluster and the mean cosmic density at the redshift of the cluster is given by \cite{kitayama+suto96}:

\begin{equation}
\Delta_{c}(z)=
18\pi^2(1+0.4093\,\omega(z)^{0.9052}) \, ,
\label{Dc}
\end{equation}
where $\omega(z)\equiv \Omega_{f}(z)^{-1}-1$ with:

\begin{equation}
\Omega_{f}(z)
=\frac{\Omega_{m,0}(1+z)^3}{\Omega_{m,0}(1+z)^3+\Omega_{\Lambda}} \, .
\label{Omz}
\end{equation}

\bigskip

\noindent The virial mass, $M_v$, and the virial radius are thus related by:

\begin{equation}
R_{v}=\left [\frac{3M_{v}}{4\pi\Delta_{c}(z)\rho_{m}(z)}
\right]^{1/3}
\label{Rv}
\end{equation}
where $\rho_{m}(z)=2.78\times10^{11}\,\Omega_{m,0}\,(1+z)^3~h^2$ M$_{\odot}$ Mpc$^{-3}$ is the mean mass density of the Universe at redshift $z$. The ratio of the gas density at high and low redshift can be assumed to scale as the ratio of the two virial densities, i.e., $\frac{\rho_{v,z=0.7}}{\rho_{v,z=0.2}}=\frac{M_{v,z=0.7}}{\left (R_{v,z=0.7} \right)^3}\times \frac {\left (R_{v,z=0.2} \right )^3}{M_{v,z=0.2}}$. In the case $M_{v,z=0.7}=M_{v,z=0.2}$, the density ratio becomes simply proportional to $\left (R_{v,z=0.2}/R_{v,z=0.7} \right )^3$. Specifically, considering both at low and high-$z$ $M_v\simeq 1.0\times 10^{15}$ M$_{\odot}$ and $\Delta M_v\simeq M_{v}/3$, the above equations give a ratio $\left (v_{i,z=0.7}/v_{i,z=0.2} \right)^3\sim 1.5$ and $\rho_{v,z=0.7}/\rho_{v,z=0.2}\sim 2.2$. The resulting ratio of turbulent energy flux derived at high and low redshift is then $\sim 3.3$.

\bigskip

\noindent The above calculation is approximate, because we only took kinematic effects and $\Lambda$CDM cosmology into account. We thus also use simulations in \cite{casano+brunetti05} to calculate the injection of turbulence during simulated cluster merging histories, and derived the turbulent energy injected in the volume swept by the infalling subcluster. We found that in this case the distributions of the turbulent energy fluxes derived for clusters with $M_v\simeq 0.9-1.4\times 10^{15}\,M_{\odot}$ in the two redshift ranges $0.2-0.3$ and $0.7-0.8$ differ by a factor of  $\sim 3-4$, in a good agreement with the simpler derivation reported above.

\bigskip
\bigskip

\noindent\textbf{Radiative lifetime and electrons acceleration in high-$z$ clusters.}
In the paper we have assumed that $\eta_{\rm rel}$, i.e. the fraction of turbulent energy flux that is absorbed by relativistic electrons (Eq. 1), is constant with redshift. This depends on the energy budget of the population of re-accelerated seed electrons, and on the microphysics of the nonlinear interaction between turbulence and particles. The latter is assumed to be redshift-independent (e.g. \cite{brunetti+lazarian07}). Thus, a constant $\eta_{\rm rel}$ implies a similar energy budget of the seed electrons that are accumulated in high- and low-$z$ ICM.

\bigskip

\noindent More specifically, we used Eq. 1 to infer a lower limit to the ratio of the magnetic field in high- and low-$z$ radio halos, thus our conclusions may change in the case that clusters at higher redshift host a population of seed electrons with a much larger energy budget. In this section we show that this possibility is very unlikely, as the higher injection rate of electrons in the ICM expected from high-redshift sources is balanced by the stronger electron energy losses.

\bigskip

\noindent The budget of seed electrons that accumulate in the ICM can be estimated as $N_e \propto Q_e \times T_{\rm max}$ (with $Q_e$ the injection rate and $T_{\rm max}$ the maximum electron lifetime). The lifetime of relativistic electrons in the ICM, which is determined by the energy losses due to Coulomb interaction with the thermal plasma and to Inverse Compton scattering with the CMB photons and synchrotron losses, can be expressed as \cite{brunetti+jones14}:

\begin{equation}
T=4\left \{\frac{1}{3}\tilde{\gamma}\beta_{z}+\frac{\tilde{n}}{\tilde{\gamma}}\left[1.2+\frac{1}{75}\ln \left (\frac{\tilde{\gamma}}{\tilde{n}}\right )\right]\right\}^{-1} \, ,
\label{eq:Tlife}
\end{equation}
where we put $\tilde{\gamma}=\gamma/300$, $\tilde{n}= n/10^{-3}$, $\beta_z=\left ( \frac{B}{3.2} \right )^2+(1+z)^4$. Here, $\gamma$ is the Lorentz factor, $n$ the number density in cm$^{-3}$ and $B$ the magnetic field strength in $\mu$G. This lifetime increases with the energy $\gamma$ of relativistic electrons at lower energy, it reaches a maximum and then decreases as a function of $\gamma$.

\bigskip

\noindent The maximum of $T$ is obtained for $\tilde{\gamma}\approx\sqrt{\frac{3\tilde{n}}{\beta_z}}$. For $B< B_{\rm CMB}$, it is determined by the combination of IC and Coulomb losses. Replacing this in Eq. \ref{eq:Tlife} gives the typical lifetime of relativistic electrons that can be accumulated in the ICM on cluster lifetime and then re-accelerated during mergers:

\begin{equation}
T_{\rm max}\,[{\rm Gyr}]=4\left \{\frac{\beta_z}{3}\sqrt{\frac{3\tilde{n}}{\beta_z}}+\tilde{n}\sqrt{\frac{\beta_z}{3\tilde{n}}}\right\}^{-1} \, .
\label{eq:Tmax}
\end{equation}
Considering that between $z=0.2$ and $z=0.7$ the virial density (and hence the thermal gas density $n$) increases by a factor of $\sim2$, the maximum lifetime ratio for electrons radiating at these two redshifts is $\sim 2.7$, considering $B \sim 1-2 ~\mu$G.

\bigskip

\noindent Several sources can inject electrons in the ICM including AGN, galaxies (galactic winds, GW) and cosmological shocks.
Most of the energy of the relativistic plasma of GW and shocks is in the form of supra-thermal and relativistic protons with only a negligible fraction in the form of electrons (e.g. \cite{brunetti+jones14}).

\bigskip

\noindent AGN are thought to be the dominant sources of relativistic electrons in the ICM (e.g. \cite{brunetti+jones14}); this hypothesis is also supported by recent LOFAR observations, where sources of fossil radio plasma from old AGN-electrons are commonly found in galaxy clusters \cite{degasperin+17,mandal+20}. The injection rate of electrons from a number of AGN in clusters, $N_{\rm AGN}$, is $Q_e \sim N_{\rm AGN} n_{\rm rel} V / \Theta$, where $n_{\rm rel}$ is the number density of radio-emitting electrons in the jets and lobes, $V$ is the volume of the radio plasma and $\Theta$ is the AGN life-time. The number density of radio-emitting  electrons in the AGN lobes can be estimated assuming equipartition between particles and magnetic fields, i.e. $n_{\rm rel}\propto \sqrt{L_{\rm syn} / V}$, with $L_{\rm syn}$ the AGN synchrotron luminosity. This gives $Q_e \propto \sqrt{L_{\rm syn} V} N_{\rm AGN}/\Theta$. The linear size of radio galaxies is observed to decline with redshift, probably due to the higher density of the medium in which radio lobes and jets expand \cite{neeser+95,blundell+99}. Thus we attempt to account for the expected change of $V$ with redshift assuming approximate pressure equilibrium of the lobes with the surrounding medium. This would imply $P_{\rm ICM} \propto n_{\rm rel} \propto \sqrt{L_{\rm syn} / V}$ and consequently $Q_e \propto L_{\rm syn} N_{\rm AGN} / ( P_{\rm ICM} \Theta )
\approx \rho_L (1+z)^3/ ( P_{\rm ICM} \Theta ) $, where $\rho_L$ is the luminosity density of radio AGN per comoving-volume. The luminosity density increases by a factor $\sim 2-2.5$ from $z\sim 0.2$ to $z\sim 0.7$ (\cite{smolciv+17} their Fig. 5). The increase of $P_{\rm ICM}$ with redshift can be estimated using virial quantities, i.e. $P_{\rm ICM} \propto n_{\rm ICM} T_{\rm ICM} \propto \frac{M_v}{R_v^3}\frac{M_v}{R_v}$, implying that pressure increases by a factor $\sim 2.5$ from $z\sim 0.2$ to $z\sim 0.7$ for clusters with the same virial mass.
We find that $Q_e$ is only about 2-3 times larger in clusters at $z\sim0.7$ if we assume that the life-time of AGN does not change much with redshift. Combining this result with the maximum lifetime of the radiating electrons, that is $\sim 2.7$ times longer at lower redshift, suggests that the budget of seed electrons at low and high redshift is similar.

\bigskip

\noindent Seed electrons in the ICM can also be generated by the decay chain of inelastic proton-proton collisions (e.g. \cite{brunetti+lazarian11,pinzke+17}). The importance of this channel depends on the energy budget of CR in the ICM, which is poorly known. A relevant contribution of re-accelerated secondary particles to radio halos is still allowed by current gamma-ray limits from FERMI-LAT (e.g. \cite{brunetti+17}). For this reason, we also compare the budget of secondary electrons accumulated in galaxy clusters at low and high redshift considering cosmological shocks as the main sources of relativistic protons. In this case, numerical simulations suggest the ratio of CR to thermal ICM pressures ($X=P_{\rm CR}/P_{\rm ICM}$) in high-$z$ clusters is generally smaller than that at low redshift (e.g. \cite{vazza+14}). The ratio of the energy budget of secondary electrons at high and low redshift can be approximately estimated as:
\begin{equation}
\frac{N_{e,z=0.7}}{N_{e,z=0.2}} \sim
\frac{X_{z=0.7}}{X_{z=0.2}} \frac{ n_{{\rm ICM},z=0.7} P_{{\rm ICM,}z=0.7}}{n_{{\rm ICM},z=0.2} P_{{\rm ICM,}z=0.2}} \frac{T_{{\rm max},z=0.7}}{T_{{\rm max},z=0.2}} \, ,
\label{secondari}
\end{equation}
implying $N_{e,z=0.7}/N_{e,z=0.2}\sim 0.4-1.8$ if we assume typical CR pressure ratios measured in simulations, i.e. $\sim 0.2-1$ (\cite{vazza+12,vazza+14}, Figs. 14 and Fig. 12 respectively), and the ratios between physical quantities as derived above: $T_{{\rm max},z=0.7}/T_{{\rm max},z=0.2}\sim 1/2.7$, $P_{{\rm ICM},z=0.7}/P_{{\rm ICM},z=0.2} \sim 2.5$ and  $n_{{\rm ICM},z=0.7}/n_{{\rm ICM},z=0.2} \sim 2$.
Similarly to the case of AGN, this implies that the budget of secondary particles that is available at higher redshift is very similar to that at low redshift.

\bigskip
\bigskip

\noindent\textbf{Reynolds number estimation.}
In Figure 3b
we have used a simple model of magnetic field amplification to infer combined constraints on the initial magnetic field and on the Reynolds number.

\bigskip

\noindent We follow the method as outlined in \cite{beresnyak12}. We assume that the amplification initially operates in a kinematic regime, where the magnetic field grows exponentially with time:

\begin{equation}
B^2(t) \sim B^2_0 \exp \left (t \, \Gamma \right)  \, ,
\label{bexp}
\end{equation}
where the time-scale of the magnetic growth is $\Gamma^{-1} \sim 30 L_d/\delta v_d$. This depends on the turbulent eddy turnover time
at the viscous dissipation scale, $L_d \sim L_{\rm inj}\, {\rm Re}^{-3/4}$,  $\delta v_d \sim v_t (L_d/L_{\rm inj})^{1/3}$. The factor $\sim 30$ is derived from simulations \cite{cho14,beresnyak+miniati16},  and accounts for the less effective stretching of the field lines due to the turbulent diffusion. When the magnetic and the kinetic energy densities become comparable at the viscous dissipation scale ($B(t)^2/8 \pi \sim 1/2 \rho \delta v_d^2$), the turbulent dynamo transits to a phase  where the magnetic field grows linearly with time. This phase is known as the non-linear regime because the magnetic field (Reynolds stress) becomes dynamically important. From Eq. \ref{bexp}, this transition occurs at a time:

\begin{equation}
t_* \sim 60 \frac{L_{\rm inj}}{v_t} {\rm Re}^{-1/2} \ln  \left ( \frac{B_*}{B_0} \right )
\label{timeexp}
\end{equation}
where $B_* \sim \sqrt{4 \pi \rho}\, v_t {\rm Re}^{-1/4}$.
According to simulations, in the non-linear regime, i.e. for time $t \geq t_*$, we assume that a constant fraction of the turbulent kinetic energy flux, $\eta_B$, is channeled into magnetic field energy ($\eta_B=0.05$ \cite{beresnyak12}). Therefore, the evolution of the magnetic field with time is:

\begin{equation}
B^2(t) \sim B^2_* + 4 \pi \frac{\rho v_t^3}{L_{\rm inj}} \eta_B (t - t_*) \, .
\label{evolution}
\end{equation}

\bigskip

\noindent In the calculations shown in Figure 3b,
we assume a continuous injection of turbulence with $v_t = 500$ km s$^{-1}$ at the injection scale $L_{\rm inj}=1$ Mpc, and $n=\rho/m_p=3\times10^{-3}$ cm$^{-3}$ (with $m_p$ the proton mass). We note that a situation of continuous injection of turbulence with these parameters is appropriate for a dynamically active and massive cluster \cite{dolag+05,donnert+18,hitomicoll18,markevitch+vikhlinin07} and consequently calculations would overestimate the magnetic field strength if extrapolated for a time scale that is much larger than that of a typical merger.

\noindent
An important parameter in the model is the time when the turbulent dynamo starts, as it fixes the time available for magnetic field amplification up to the epoch of observation. Here, we assume that the turbulent dynamo starts when the clusters have formed a quarter of their mass. This implies that the time available for the amplification is of about 3.7~Gyr \cite{fakhouri+10,giocoli+12}. Assuming a turbulent velocity of $v_t = 500$~km~s$^{-1}$ at an earlier stage in the cluster life-time is optimistic, as the cluster would be too small. Since the turbulent energy flux is proportional to $v_t^3$, the effect of much weaker turbulence at earlier epochs is not relevant for the magnetic field evolution. Conversely, if we assume that the dynamo started later in the cluster life-time, e.g. when it has assembled half of its mass, it results in more stringent Reynolds numbers, since the time available for the  amplification is shorter (about 2.7~Gyr). Lower turbulent velocities and lower number densities would imply higher values of the Reynolds number. On the other hand, higher turbulent velocities would imply a smaller value of Re. However, even if we consider $v_t=800$ km s$^{-1}$ \cite{dolag+05,donnert+18}, i.e. the case of a large turbulent pressure of about 30--40\% the ICM pressure, we obtain Reynolds numbers that are only three times smaller than the case with $v_t = 500$ km s$^{-1}$. We also mention that in real clusters turbulence is induced in a medium that is highly stratified due to gravity \cite{roh+19}. In this case, the transport of magnetised turbulent eddies toward the outskirts could make the turbulent dynamo slightly less efficient than in our model, implying that our constraints on the Reynolds number are conservative.