\section{Additional Details of The Proposed Method}

Figure~\ref{fig:rotation_acgan} depicts the proposed method. Real and fake/generated images are all rotated four directions. For the shared discriminator D, its goal is to predict true/fake labels of the \textit{non-rotated} images, and detect the rotation of the \textit{real} images.
For the generator, its goal is to generate images to fool discriminator, at the same time trying to generate images that are easier for discriminator to detect rotations (adjusted by $\alpha$).

\begin{figure}[h]
\includegraphics[width=\textwidth]{figure/ssgan_model.pdf}
\caption{\label{fig:rotation_acgan} Rotation feature learning discriminator. The dotted line arrow uses only the non-rotated images, for real/fake classification task. The solid line arrow uses all the four rotated images, for the rotation classification task. }
\end{figure}

\section{Further Experiments and Details}

\paragraph{Datasets}
We consider four datasets: \textsc{cifar10}{}, \textsc{celeba-hq}{}, \textsc{lsun-bedroom}{}, and
\textsc{imagenet}{}.  The \textsc{lsun-bedroom}{} dataset~\citep{yu15lsun} contains around 3M
images.  We partition the images randomly into a test set containing
30588 images and a train set containing the rest.  \textsc{celeba-hq}{}
contains 30k images~\citep{karras2017progressive}.  We use the
$128\times128\times3$ version obtained by running the code provided by
the authors\footnote{Available at
\url{https://github.com/tkarras/progressive\_growing\_of\_gans}.}.  We
use 3000 examples as the test set and the remaining examples as the
training set.  \textsc{cifar10}{} contains 70K images ($32\times32\times3$),
partitioned into 60000 training instances and 10000 testing instances.
Finally, we evaluate our method on \textsc{imagenet}{}, which contains
$1.3$M training images and $50$K test images.  We re-size the images
to $128\times128\times3$ as done in~\citet{miyato2018cgans}
and~\citet{zhang2018self}.

\paragraph{Metrics}
To quantitatively evaluate generated samples from different methods, we mainly \emph{Frechet Inception Distance} (FID)~\citep{heusel2017gans}. In FID, the true data and generated samples are first embedded in some feature space (a specific layer of Inception).
Then, a multivariate Gaussian is fit to the data and the distance computed as
$\FID(x, g) = ||\mu_x -\mu_g||_2^2 + \Tr(\Sigma_x + \Sigma_g - 2(\Sigma_x\Sigma_g)^\frac12)$,
where
$\mu$ and $\Sigma$ denote the empirical mean and covariance
and subscripts $x$ and $g$ denote the true and generated data respectively.
As shown in \citet{sajjadi2018assessing,lucic2018}, FID is shown to be sensitive to both the addition of spurious modes and to mode dropping.

\paragraph{Additional experimental details}
When the rotation loss is added, we take a quarter of images in the batch, rotate for all four directions to synthetic images of known rotation orientations. This subset rotation trick could reduce and keep the extra computation small. It is worth noting that other than the extra loss term, other parts (architectures, hyper-parameters) are all kept the same for our method and baselines (both unconditional and conditional GANs). For simplicity, we fix $\beta=1$ and compare $\alpha$ in $\{0.2, 0.5, 1\}$; $\alpha=0.2$ is used by default.

We use the Adam optimizer. The learning rate is fixed for both datastes as $0.0002$. All models are trained for $100$k steps\footnote{Results in the main paper for \textsc{imagenet}{} are obtained from training with 1M steps.}. All models are trained with batch size of 64 in a single Nvidia P100 GPU. Reported results by default are the mean of from 3 independent runs with different random seeds.

\paragraph{Robustness test
We test the effects of training stabilization brought by the SsGAN.
We consider two types of hyper-parameter settings: First, controlling the Lipschitz constant of the discriminator, which is a central quantity analyzed in the GAN literature~\citep{ miyato2018spectral,zhou2018understanding}.
We
consider two state-of-the-art techniques: Gradient Penalty~\citep{gulrajani2017improved}, and Spectral Normalization~\citep{miyato2018spectral}.
For the Gradient Penalty regularizer we consider regularization strength $\lambda \in \{1,10\}$.
Second, we consider different Adam optimizer hyper-parameters.
We test two popular settings $(\beta_1, \beta_2)$: $(0.5, 0.999)$ and $(0,0.9)$.  Previous studies find that multiple discriminator steps per generator step help training~\citep{goodfellow2014generative,salimans2016improved}, thus we also consider both $1$ and $2$ discriminator steps per generator step\footnote{We also experimented with 5 steps which didn't outperform the $2$ step setting.}.  In total, this amounts to three different sets of hyper-parameters for $(\beta_1, \beta_2,
\text{disc\_iter})$: $(0, 0.9, 1)$, $(0, 0.9, 2)$, $(0.5, 0.999, 1)$.

\subsection{Quantitative Comparisons of the Generative Samples}
\label{fid_quantitative_compare}

\paragraph{Convergence curve} Figure \ref{fig:convergence_cruves} shows the convergence curves of different methods. We can clearly see the proposed SsGAN significantly improves the performance over the other unconditional GAN, SN-GAN, and even matches the conditional P-cGAN.

\begin{figure}[ht]
  \centering
  \begin{subfigure}[b]{0.425\textwidth}
    \includegraphics[width=1\textwidth]{figure/convergence_fid_cifar10.pdf}
    \caption{\textsc{cifar10}{}}
  \end{subfigure}
  \begin{subfigure}[b]{0.4\textwidth}
    \includegraphics[width=1\textwidth]{figure/convergence_fid_imagenet_1m.pdf}
    \caption{\textsc{imagenet}{}}
  \end{subfigure}\hfill
  \caption{\label{fig:convergence_cruves} FID convergence curves for different methods with 100k training steps. The self supervision close the gap between the conditional and unconditional GANs, sometimes even surpass the class label conditional GAN.}
\end{figure}

\paragraph{Robustness improvement} The use of self-supervised loss not only improves the performance of the default hyper-parameters we use, but also improves under multiple types of hyper-parameters. Table \ref{tab:gp_robustness_fid} shows the performance comparison with Gradient Penalty with various hyper-parameters, and Table \ref{tab:sn_robustness_fid} shows the performance comparison with Spectral Normalization with various hyper-parameters. Both figures demonstrate SsGAN reduces the sensitivity of GAN training to hyper-parameters, hence the training is more robust.

\begin{table}
\centering
\small
\caption{\label{tab:gp_robustness_fid} FID comparisons for unconditional GANs under various hyper-parameter settings with gradient penalty. Self-supervision loss reduces the sensitivity of GAN training to hyper-parameters.}
\begin{tabular}{cccccccc}
\toprule
   &     &       &  & \multicolumn{2}{c}{\textsc{cifar10}{}} & \multicolumn{2}{c}{\textsc{imagenet}{}} \\
$\lambda$ & $\beta_1$ & $\beta_2$ & Disc\_iters  &      SN-GAN &   SsGAN &      SN-GAN &    SsGAN \\
\midrule
\multirow{3}{*}{1} & \multirow{2}{*}{0.0} & \multirow{2}{*}{0.900} & 1 &  121.047 &  \textbf{25.798} &  183.364 &   \textbf{80.669} \\ \cline{4-8}
   &     &       & 2 &   28.107 &  \textbf{26.984} &   85.130 &   \textbf{83.081} \\
\cline{2-8}
\cline{3-8}
   & 0.5 & 0.999 & 1 &   78.542 &  \textbf{25.890} &  104.731 &   \textbf{91.628} \\
\cline{1-8}
\multirow{3}{*}{10} & \multirow{2}{*}{0.0} & \multirow{2}{*}{0.900} & 1 &  188.523 &  \textbf{28.481} &  227.041 &   \textbf{85.376} \\ \cline{4-8}
   &     &       & 2 &   29.110 &  \textbf{27.735} &  227.738 &   \textbf{80.822} \\
\cline{2-8}
\cline{3-8}
   & 0.5 & 0.999 & 1 &  117.671 &  \textbf{25.216} &  242.705 &  \textbf{144.353} \\
\bottomrule
\end{tabular}
\end{table}

\begin{table}
\centering
\small
\caption{\label{tab:sn_robustness_fid} FID comparisons for unconditional GANs under various hyper-parameter settings with spectral normalization. Self-supervision loss reduces the sensitivity of GAN training to hyper-parameters.}
\begin{tabular}{cccccccc}
\toprule
    &       &  & \multicolumn{2}{c}{\textsc{cifar10}{}} & \multicolumn{2}{c}{\textsc{imagenet}{}} \\
$\beta_1$ & $\beta_2$ & Disc\_iters  &     SN-GAN &   SsGAN &      SN-GAN &   SsGAN \\
\midrule
\multirow{2}{*}{0.0} & \multirow{2}{*}{0.900} & 1 &  87.864 &  \textbf{19.645} &  129.964 &  \textbf{86.095} \\ \cline{3-7}
    &       & 2 &  20.245 &  \textbf{17.885} &   80.047 &  \textbf{70.637} \\
\cline{1-7}
\cline{2-7}
0.5 & 0.999 & 1 &  86.869 &  \textbf{18.235} &  201.936 &  \textbf{99.968} \\
\bottomrule
\end{tabular}
\end{table}

\paragraph{Selection of the hyper-parameters} We further compare different choices of $\alpha$ (while fixing $\beta=1$). As discussed, an inappropriate value of $\alpha$ could bias the generator, while a reasonable one would help. As shown in Table \ref{fig:beta_variant}, we indeed see the effectiveness of $\alpha$. In the values compared, the optimal $\alpha$ is 1 for \textsc{cifar10}{}, and 0.2 for \textsc{imagenet}{}. By default we use 0.2 for all datasets for simplicity.

\begin{figure}[h]
  \centering
  \begin{subfigure}[b]{0.26\textwidth}
    \includegraphics[width=1\textwidth]{figure/alpha_variant_fid_score_cifar10.pdf}
    \caption{\textsc{cifar10}{}}
  \end{subfigure}
  \begin{subfigure}[b]{0.25\textwidth}
    \includegraphics[width=1\textwidth]{figure/alpha_variant_fid_score_imagenet.pdf}
    \caption{\textsc{imagenet}{}}
  \end{subfigure}
  \caption{\label{fig:beta_variant} Performance under different $\alpha$ values.}
\end{figure}

\subsection{Qualitative Comparisons of the Generative Samples}

Figure \ref{fig:samples_imagenet} shows the generated examples from different GANs for \textsc{imagenet}. Although the baseline unconditional GAN generates sharp images, majority of the generated objects e.g. dogs are distorted. With the supervised labels, P-cGAN's samples are improved, e.g. we see the dog faces and diverse objects. SsGAN generates similar quality images to the conditional P-cGAN.

\begin{figure}[ht]

\begin{subfigure}[b]{0.49\textwidth}
    \begin{center}
     \includegraphics[width=1\textwidth]{figure/imgs_imagenet_1m_sngan}
    \end{center}
    \caption{SN-GAN (FID=59.3)}
\end{subfigure}\hfill
\begin{subfigure}[b]{0.49\textwidth}
    \begin{center}
     \includegraphics[width=1\textwidth]{figure/imgs_imagenet_1m_pcgan}
    \end{center}
    \caption{P-cGAN (FID=43.0)}
\end{subfigure}\hfill
\begin{subfigure}[b]{0.49\textwidth}
    \begin{center}
     \includegraphics[width=1\textwidth]{figure/imgs_imagenet_1m_ssgan}
    \end{center}
    \caption{SsGAN (FID=47.8)}
\end{subfigure}\hfill
\begin{subfigure}[b]{0.49\textwidth}
    \begin{center}
     \includegraphics[width=1\textwidth]{figure/imgs_imagenet_1m_ssgan_sbn}
    \end{center}
    \caption{SsGAN (sBN) (FID=43.9)}
\end{subfigure}\hfill
\caption{\label{fig:samples_imagenet} Random samples generated by GANs described in this paper on \textsc{imagenet} (at 1M steps).}
\end{figure}

\subsection{Comparison of Representation Quality}

We compare the representation qualities of the discriminator's intermediate layers. The comparison is performed by training logistic regression classifiers on top of the feature maps from each ResNet block to perform the 1000-way \textsc{imagenet}{} classification task, as proposed in \citet{zhang2016colorful}. The representation features are spatially resized
to around 9000 dimension with adaptive max pooling. We report top-1 classification accuracy.

We compare the SN-GAN, P-cGAN and SsGAN (SBN) models with the best FID score from Section \ref{fid_quantitative_compare}. We also consider ablating the GAN loss from SsGAN, to yield a rotation only model (Rotation). Rotation uses the same architecture and hyper parameters as the SsGAN discriminator.
Evaluation is performed on three independent models with different random seeds, where the mean accuracy and the standard deviation are reported. We tune the learning rate for linear regression on a validation set created by holding out $10\%$ of the training set. The linear regression results are insensitive to the hyperparameters. On \textsc{cifar10}{}, we use SGD with batch size 256, momentum 0.9 and learning rate 0.1. We drop the learning rates by a factor of 10 after every 30 epochs. We train in total for 100 epochs, where the learning rates are dropped by a factor of 10 after every 30 epochs.
On \textsc{imagenet}{}, we use SGD with batch size 64, momentum 0.9 and learning rate 0.01. We train in total for 30 epochs, where the learning rates are dropped by a factor of 10 after every 10 epochs.

Table~\ref{tab:cifar10_top1} shows the top-1 classification accuracy on \textsc{cifar10}{}. Figure~\ref{fig:accuracy-steps-cifar} shows the results as s function of training steps of the original model. We plot the quality of representations extracted from each of the the 4 blocks for ResNet \textsc{cifar10} architecture. The conditional model, P-cGAN, produces  similar results to the other unsupervised model GAN and Rotation model. The SsGAN outperforms the other models under most conditions.

Table \ref{tab:imagenet_top1} and Figure~\ref{fig:accuracy-steps-imagenet} show results on \textsc{imagenet}{}.
There are six blocks in the \textsc{imagenet}{} Resnet architecture.
For SN-GAN and Rotation, block3 performs best, for P-cGAN and SsGAN, block 5 performs best.
It is expected that P-cGAN would benefit from representations closer to the classification layer because it sees the actual labels being used in the downstream representation evaluation.
In block5, we observe that SsGAN and P-cGAN representation quality is not improving after 200K steps, while SsGAN improves the representation quality for all 1M training steps.
Overall, representations are improved with self-supervised GAN and conditional GAN, which correlated with their improved samples.

Table~\ref{tab:imagenet_compare} shows the comparison of our method to other self supervised learning methods. We select the best run according to the validation results. Overall, the SsGAN model achieves comparable results on \textsc{imagenet}{}. Note that our method is not tuned for representation, but rather for image generation. Among those methods, BiGAN~\citep{donahue2016adversarial} learns the representation from GANs as well. BiGAN learns the representation with an additional encoder network, while our method is arguably simpler because it extracts the representation directly from the discriminator network.

\begin{table}[h]
\centering
\caption{\label{tab:cifar10_top1} Top-1 accuracy on \textsc{cifar10}.}

\begin{tabular}{llllll}
\toprule
Method & SN-GAN & P-cGAN & Rotation & SsGAN (sBN) \\
\midrule
Block0    &   $0.719 \pm 0.002$ &   $0.719 \pm 0.003$ & $0.710 \pm 0.002$ & $0.721 \pm 0.002$\\
Block1      &   $0.762 \pm  0.001$ &   $0.759 \pm  0.003$ & $0.749 \pm 0.003$ & $0.774 \pm 0.003$ \\
Block2 &   $\mathbf{0.778 \pm  0.001}$ &   $0.776 \pm  0.005$ & $\mathbf{0.762 \pm 0.003}$ & $0.796 \pm 0.005$ \\
Block3  &  $0.776 \pm  0.005$ &   $\mathbf{0.780 \pm  0.006}$ & $0.752 \pm 0.006$ & $\mathbf{0.799 \pm 0.003}$ \\
\bottomrule
\end{tabular}
\end{table}

\begin{table}[h]
\centering
\small
\caption{\label{tab:imagenet_top1} Top-1 accuracy on \textsc{imagenet}{} for all Resnet blocks.}

\begin{tabular}{lllll}
\toprule
Method & SN-GAN & P-cGAN & Rotation & SsGAN (sBN) \\
\midrule
Block0      &    $0.084 \pm 0.076$   &    $0.169 \pm  0.001$    &    $0.172 \pm  0.002$    &   $0.170 \pm  0.001$ \\
Block1      &    $0.064 \pm  0.107$   &    $0.190 \pm  0.004$    &    $0.169 \pm  0.005$    &   $0.198 \pm  0.003$ \\
Block2      &    $0.078 \pm  0.134$   &    $0.232 \pm 0.002 $    &    $0.202 \pm 0.004$     &   $0.250 \pm 0.001$ \\
Block3      &    $\mathbf{0.087 \pm 0.149} $   &    $0.271 \pm 0.003$    &    $\mathbf{0.239 \pm 0.005}$     &   $0.284 \pm 0.002$ \\
Block4      &    $0.0744 \pm 0.127$   &    $0.273 \pm 0.019  $    &    $0.180 \pm 0.003$     &   $0.345 \pm 0.003$ \\
Block5      &    $0.075 \pm 0.128$   &    $\mathbf{0.324 \pm 0.007}$    &    $0.176 \pm 0.033$     &   $\mathbf{0.353 \pm 0.006}$ \\
\bottomrule
\end{tabular}
\end{table}

\begin{table}[h]
\centering
\small
\caption{\label{tab:imagenet_compare} Comparison with other self-supervised representation learning methods by top-1 accuracy on \textsc{imagenet}{}.}

\begin{tabular}{llllll}
\toprule
Method & Context~\citep{doersch2015unsupervised} & BiGAN~\citep{donahue2016adversarial} & Colorization~\citep{zhang2016colorful} & RotNet~\citep{gidaris2018unsupervised} & SsGAN (sBN) \\
\midrule
Accuracy      &    0.317   &    0.310    &    0.315    &   0.387 & 0.359 \\
\bottomrule
\end{tabular}
\end{table}

\begin{figure}[h]
\begin{center}
\epsfig{file=figure/cifar_step_to_accuracy,height=8.0cm}
\end{center}
\caption{\label{fig:accuracy-steps-cifar} Top 1 accuracy on \textsc{cifar10}{} with training steps from 10k to 100k.}
\end{figure}

\begin{figure}[h]
\begin{center}
\epsfig{file=figure/imagenet_step_to_top1_accuracy,height=12.0cm}
\end{center}
\caption{\label{fig:accuracy-steps-imagenet} Top 1 accuracy on \textsc{imagenet}{} validation set with training steps from 10k to 1M.
}
\end{figure}

\section{Experiments}
\label{sec:experiments}

We show that self-supervision improves GAN training, even matching the performance of conditional GANs.
We also demonstrate, by way of representation evaluation, that the features learned by the self-supervised GAN are more meaningful than those of a vanilla GAN.

\paragraph{Datasets and Settings}
Our main result uses the \textsc{imagenet} dataset.
For the self-supervised GAN (SsGAN), we adopt the architectures of, and compare to, the unconditional GAN in~\citet{miyato2018spectral} (SN-GAN) and conditional GAN in~\citet{miyato2018cgans} (P-cGAN).
Both the generator and discriminator consist of ResNet architectures.
For the conditional generator in P-cGAN, we apply label-conditional batch norm.
SsGAN does not use labels, so by default the generator does not use condition batch-norm, as in SN-GAN.
Therefore, we also consider a variant of the SsGAN with self-modulated batch norm (sBN)~\citep{abn} on the generator to replace label-conditional batch norm.
We compare image generation quality using the Frechet Inception Distance (FID)~\citep{heusel2017gans}.
We train the model for 1M steps using a single P100 GPU.
Further details and results can be found in the supplementary material.

We also test three smaller datasets, \textsc{cifar10}, \textsc{lsun-bedroon}, and \textsc{celeba-hq}.
Labels for Pc-GAN are only available for \textsc{cifar10}.
For all datasets we set $\beta=1$ and
$\alpha=0.2$, chosen using a small hyperparameter sweep.
For all other hyperparameters, we use the values in \citet{miyato2018spectral} and \citet{miyato2018cgans}.
We train on these datasets for $100$k steps on a single P100 GPU.

\paragraph{Image generation}
Figure~\ref{fig:convergence_fid_imagenet_1m} shows FID training curves on
\textsc{imagenet}.
The unconditional GAN is unstable, and sometimes diverges.
The conditional P-cGAN substantially outperforms the unconditional SN-GAN.
The SsGAN is stable, and even performs as well as the P-cGAN.

\begin{minipage}{\textwidth}
\centering
\begin{minipage}[b]{0.53\textwidth}
    \centering
    \includegraphics[width=\linewidth]{figure/convergence_fid_imagenet_1m}
    \captionof{figure}{\label{fig:convergence_fid_imagenet_1m}Mean training curves across three random seeds.}
\end{minipage}
\hfill
\begin{minipage}[b]{0.42\textwidth}
\vspace{0pt}
    \centering
    \small
    \begin{tabular}{lll}
    \toprule
    Dataset & Method &   FID   \\
    \midrule
    \multirow{4}{*}{\textsc{cifar10}} & SN-GAN &     19.725 \\
            & P-cGAN &     \textbf{15.604} \\
            & SsGAN &     17.114 \\
            & SsGAN (sBN) &     \textbf{15.653} \\ \midrule
    \multirow{4}{*}{IMAGENET} & SN-GAN &     56.670 \\
            & P-cGAN &     \textbf{42.071} \\
            & SsGAN &     47.555 \\
            & SsGAN (sBN) &     43.869 \\ \midrule
    \multirow{3}{*}{BEDROOM} & SN-GAN &     16.021 \\
            & SsGAN &     13.664 \\
            & SsGAN (sBN) &     \textbf{13.298} \\ \midrule
    \multirow{3}{*}{CELEBA-HQ} & SN-GAN &     \textbf{23.765} \\
            & SsGAN &     26.113 \\
            & SsGAN (sBN) &     24.355 \\
    \bottomrule
    \end{tabular}
    \captionof{table}{\label{tab:best_fid} Best FID attained across three random seeds.}
\end{minipage}
\end{minipage}

\begin{figure}[ht]
  \centering
  \begin{subfigure}[b]{0.48\textwidth}
    \includegraphics[width=1.0\textwidth]{figure/cifar_step_to_accuracy_block3}
    \caption{\label{fig:accuracy-steps-cifar-block3}CIFAR-10}
  \end{subfigure}
  \begin{subfigure}[b]{0.48\textwidth}
    \includegraphics[width=1.0\textwidth]{figure/imagenet_step_to_top1_accuracy_block5.pdf}
    \caption{\label{fig:accuracy-steps-imagenet-block5}IMAGENET}
  \end{subfigure}\hfill
  \caption{\label{fig:accuracy-steps}
  Top 1 accuracy (mean across three seeds) to predict labels from discriminator representations from the final ResNet block.
  x-axis gives the number of GAN training iterations.}
\end{figure}

Table~\ref{tab:best_fid} shows FID of the best run across three random seeds for each dataset/model combination.
Again, for \textsc{cifar} the SsGAN (with sBN) improves unconditional performance, and even equals P-cGAN.
For \textsc{bedroom}, we also see some improvement.
Self supervision appears not to help in \textsc{celeba-hq},
we believe that this data is too simple and face orientations do not provide useful signal.

\paragraph{Representation Quality}
We compare each model's discriminator representations.
We also ablate the GAN loss from SsGAN to assess the representation learned from self-supervision alone.
We followed the evaluation protocol of \cite{zhang2016colorful}.
We train a logistic regression to predict the class labels from the feature maps of the discriminator
and report top-1 classification accuracy.

Figure~\ref{fig:accuracy-steps} shows the classification performance using the final ResNet block as a function of GAN training iterations.
SsGAN produces better representations than SN-GAN.
The same is true for all the other blocks (see supplementary, Figures~\ref{fig:accuracy-steps-cifar} and~\ref{fig:accuracy-steps-imagenet}).
Interestingly, during training the SsGAN overtakes P-cGAN, which sees the class labels. This indicates that P-cGAN is overfitting the training data.
Using rotation loss alone substantially decreases representation quality of the SsGAN.
Note that our model uses hyperparameters designed for image generation;
therefore, we do not seek state-of-the-art representations.
Instead, we use representation evaluation to understand and motivate our strategy.
Indeed it seems that self-supervision improves both discriminator representations and image generation.

\section{The Self-supervised GAN}

We first discuss discriminator forgetting and motivate representation learning for GANs. We then present our solution using rotation based self-supervision.

\subsection{Discriminator Forgetting}

The original value function for GAN training is:
\begin{equation}
\label{eq:gan_obj_classic}
V(G, D) = \mathbb{E}_{\bm x\sim P_{\mathrm{data}}(\bm x)}[\log P_D(S=\mathrm{real}|\bm x)]
+ \mathbb{E}_{\bm x\sim P_{G}(\bm x)}[\log (1- P_D(S=\mathrm{real}|\bm x))]
\end{equation}
where $P_{\mathrm{data}}$ is the true data distribution, and $P_G$ is the data distribution
induced by feeding noise drawn from a simple distribution $\bm z\sim P(\bm z)$ through the generator, $\bm x = G(\bm z)$.
$P_D(S|\bm x)$ is the discriminator's Bernoulli distribution over the sources (real or fake).

The generator maximizes Equation~\ref{eq:gan_obj_classic}, while the discriminator minimizes it.
Training is typically performed via alternating stochastic gradient descent.
The discriminator, therefore, classifies between two data distributions, $P_{\mathrm{data}}(\bm x)$, and $P_G^{(t)}(\bm x)$, where $t$ indexes training iterations.
The latter, however, is non-stationary because the generator's parameters are updated over time.
This is an online learning problem for which explicit temporal dependencies have been proposed to improve training~\citep{anonymous2019generative,shrivastava2017learning,grnarova2017online,salimans2016improved}.

In online learning, neural networks tend to forget previous tasks~\cite{kirkpatrick2017overcoming,mccloskey1989catastrophic,french1999catastrophic}.
We demonstrate forgetting in a toy setting, Figure~\ref{fig:non-stationary}(a).
Here, a classifier to performs sequential 1-vs.-all classification task on each of the ten classes in \textsc{cifar10}.
The classifier is trained for $1$k iterations on each task before switching to the next.
Figure~\ref{fig:non-stationary}(a) shows substantial forgetting, despite the tasks being similar.
The classifier returns almost to random accuracy, $0.5$, each time the task switches.
The accuracy even drops substantially when returning to the original task at iteration $10$k.
This demonstrates that the model does not retain generalizable representations during online learning.

Consideration of convergence provides further indication that discriminators may forget.
\citet{goodfellow2014generative} show that the optimal discriminator estimates the likelihood ratio between the generated and real data distributions.
Therefore, given a perfect generator, where $P_G = P_D$, the optimal discriminator simply outputs $0.5$, and has no requirement to retain any meaningful representations.

Discriminator forgetting may cause training difficulties because
(1) it does not learn meaningful representations to guide the generator, and
(2) the generator can revert to generating old images to fool it~\citep{thanh2018catastrophic}.
This is particularly problematic if the data is diverse and complex and the generator exhibits mode collapse.
Therefore, we add self-supervision to encourage the discriminator to retain useful representations.

\begin{figure}[ht]
\begin{center}
\epsfig{file=figure/nonstationary,height=4.5cm}
\end{center}
\caption{\label{fig:non-stationary}
Image classification accuracy when the underlying class distribution shifts every $1$k iterations.
The vertical dashed line indicates the end of an entire cycle, and return to the original classification task at $t=0$.
\emph{Left}: vanilla classifier.
\emph{Right}: classifier with additional self-supervised loss.
This example demonstrates that a classifier/discriminator may fail to learn generalizable representations in a non-stationary environment, but an auxiliary loss mitigates this problem.}
\end{figure}

\subsection{The Self-Supervised GAN}

Self-supervised learning is a family of methods for building representations from unsupervised data.
Self-supervision works by creating artificial supervised tasks from unsupervised data, training on these tasks, and then extracting representations from the resulting networks~\citep{dosovitskiy2014,doersch2015unsupervised,zhang2016colorful}.
Here, we apply the successful image rotation self-supervision method~\citep{gidaris2018unsupervised}.
In this method, the self-supervised network predicts the angle of rotation of an image.
In Figure~\ref{fig:non-stationary}(b) we motivate this loss using our toy problem.
When we add the self-supervised loss, the network learns features that transfer across tasks; performance continually improves, and does not drop to $0.5$ when the distribution shifts.

For the self-supervised GAN, the specific losses we use for the generator and discriminator are:
\begin{align}
L_G &= -V(G, D) -\alpha \mathbb{E}_{\bm x\sim P_G}[\sum_{r\in \mathcal{R}}\nolimits \log P^\prime_{D}(R=r|{\bm x}^r)],\\
L_D &= V(G, D) -\beta \mathbb{E}_{\bm x \sim P_{\mathrm{data}}}[\sum_{r\in \mathcal{R}}\nolimits \log P^\prime_{D}(R=r|{\bm x}^r)],
\end{align}
where $V(G, D)$ is the original GAN loss in Equation~\ref{eq:gan_obj_classic}.
$r\in\mathcal{R}$ is a rotation selected from a set of possible rotations.
We use
$\{ 0^{\degree}, 90^{\degree}, 180^{\degree}, 270^{\degree} \}$
as in~\citet{gidaris2018unsupervised}.
$P'_D(R|\bm x^r)$ is the discriminator's distribution over rotations,
and $\bm x^r$ is the image $\bm x$ transformed by rotation $r$.
Architecturally, we use a single discriminator network with two heads to
compute $P_D$ and $P^\prime_{D}$. See supplementary Figure~\ref{fig:rotation_acgan} for an overview.

\paragraph{A note on convergence}
With $\alpha>0$ convergence, even under optimal conditions, to the true data distribution $P_G=P_{\mathrm{data}}$ is not guaranteed.
This may not be a concern because current GANs are far from attaining the optimal solution.
If it is, one could anneal $\alpha$ to zero during training.
Our intuition is that the proposed loss encourages the discriminator to learn and retain meaningful representations that allow it to distinguish rotations as well as true/fake images.
The generator is then trained to match distributions in this feature space which encourages the generation of realistic objects.

\section{Introduction}

Generative Adversarial Networks (GANs) are a class of unsupervised generative models~\citep{goodfellow2014generative}.
GANs involve training a \emph{generator} and \emph{discriminator} network in an adversarial game, such that the generator learns to produce samples from a desired data distribution.
Training GANs is challenging because one searches for a Nash equilibrium of a non-convex game in a high-dimensional parameter space.
GANs are typically trained with alternating stochastic gradient descent.
However, this training procedure is unstable and lacks guarantees~\cite{salimans2016improved}.
Training may exhibit instability, divergence, cyclic behaviour, or mode collapse~\citep{mescheder2018training}.
As a result, many works propose techniques to stabilize GAN training~\citep{mao2016least,gulrajani2017improved,miyato2018spectral,abn,radford2016,zhang2018self,karras2017progressive}.

Recently, several papers frame GAN training as a continual, or online, learning problem~\citep{thanh2018catastrophic,anonymous2019generative,seff2017continual,shrivastava2017learning,kim2018memorization,grnarova2017online}.
GANs may be viewed as online learning because each network, the generator and discriminator, must learn in the context of a non-stationary environment.
The discriminator, for example, is a binary classifier where one class (the fake samples) is non-stationary.
With this view, \citet{thanh2018catastrophic} study catastrophic forgetting in the discriminator and mode collapse, relating these to training instability.
\citet{anonymous2019generative} counter discriminator forgetting by leveraging techniques from continual learning directly (Elastic Weight Sharing~\cite{kirkpatrick2017overcoming} and Intelligent Synapses~\citep{zenke2017continual}).
Other works address forgetting in GANs by retraining on old data, or adding an explicit memory~\citep{shrivastava2017learning,kim2018memorization}.
Finally, \citet{seff2017continual} extend continual learning for GANs by training on nonstationary true-data distributions also.

We follow this line of work, but propose an alternative approach.
Instead of adding explicit replay or memorization strategies, we apply representation learning.
We train the discriminator to learn stable representations that are useful for identifying real images.
We propose to add self-supervision to the discriminator,
in particular, we apply a rotation-based loss~\citep{gidaris2018unsupervised}.

Conditional GANs, which use side information such as class labels, are state-of-the-art in generating high fidelity complex images~\cite{mirza2014conditional,odena2017,miyato2018cgans}.
A major contributor to their success is augmentation of the discriminator with supervised information, which encourages it to learn stable and diverse representations.
However, labels are not always available, and unconditional GANs perform much worse.
Our self-supervised GAN requires no additional information, and closes the performance gap for image generation between conditional and unconditional GANs.
As intended, we observe that the self-supervised GAN learns better discriminator representations than other GANs for the task of image classification.
We hope that this work encourages further investigation into representation and continual learning to improve GANs.

\section{Related Work}

\section{Conclusion}
Based on the intuition discriminator can have a hard time learning representation in non-stationary setting, we proposed to adopt self-supervised method to enhance GAN training, which closes the performance gap between conditional and unconditional GANs.


\section{Introduction}

A vital computational problem within probabilistic graphical models such as Bayesian networks is the problem of finding the {\em mode} or {\em most probable explanation} of a set of variables given observed values for other variables in the network. When the network includes latent or hidden variables (i.e., variables that have neither been observed nor are of interest for the explanation) this problem is known as \textsc{Partial MAP}; the explanation sought is the MAP explanation, i.e., the joint value assignment to the explanation variables that has maximum posterior probability. In this paper we are interested in the problem of finding not specifically the MAP explanation, but {\em other} explanations that have desirable properties. There are two distinct reasons why we may be interested in alternative explanations:

\begin{enumerate}
\item As a means of approximating the MAP explanation. \textsc{Partial MAP}\ is a highly intractable problem \cite{Campos11,Kwisthout11} and we may find an acceptable approximation thereof (to be further explicated later) `good enough';
\item To obtain `alternative good' explanations {\em in addition to} the MAP explanation, to allow us to explore several alternatives (e.g. search results) or to cover for a set of likely explanations (e.g. medical conditions).
\end{enumerate}

`Good enough' and `alternative good' explanations are conceptually different when we look at the {\em structure} of the explanation as compared to the MAP explanation; that is, how similar or dissimilar the joint value assignments (of the MAP explanation and the alternative explanation) are. In the first case we are more than happy to obtain an explanation that is almost identical to the MAP explanation; in fact, in some problem domains this might even be a prerequisite. For example, in computational cognitive modeling an explanation that does not resemble the MAP explanation would not be acceptable as a valid approximation of the MAP explanation, even if it has a comparable probability. In the second case, our goal is to find alternative explanations that are {\em structurally different} yet are still plausible. For example, in response to a search query we don't want to end up with ten results that refer to minor variations of essentially the same web-page, even if they happen to be the most probable given the query.

Note that {\em structure} approximation (where we seek an explanation with Hamming distance at most $d$ of the MAP explanation) is really different from {\em value} approximation (where we seek an explanation with almost-as-high probability, e.g. within ratio $r$ of the MAP explanation) and from {\em rank} approximation (where we seek an explanation that is within the $m$ best explanations) \cite{Kwisthout15}. These notions are really orthogonal, as Figure \ref{map_structs} will reveal.

\makefig{map_structs}
{Graphical depiction of possible relationships between similarity and probability of explanations, with on the X-axis defines an order of the explanations according to their similarity to the best explanation. In ``linear decrease'' the structure and the probability correlate almost completely; in ``random noise'' almost not at all. There can be a sole peak of probability mass (``lonely mountain'') with all other explanations having almost zero probability. Alternatively, there can be several structurally distinct ``local maxima'' of probability mass. Particularly in the latter case it might be interesting to look at dissimilar explanations that have a relatively high probability.}
{12cm}{jpg}

Current algorithms that seek to find alternative explanations by exploring local maxima (e.g. \cite{Chen13,Felzenszwalb10,Chen18}) may fare well if the explanatory landscape is as in the ``local maxima'' panel of Figure \ref{map_structs}. They might not find good explanations, or take a lot of time, if the landscape is in either of the other panels of Figure \ref{map_structs}. However, the computational complexity of this problem has not yet been investigated. In this paper we are specifically interested in explanations that rank well (are within the best $m$ explanations) {\em and} that are either structurally similar or dissimilar; that is, we complement the results of \cite{Kwisthout15}. In this paper we will further explicate both problems and explore the computational complexity of both of them. We will start with offering some necessary preliminaries and sharing our notational conventions. In Section \ref{Main_results} we will formalize both problems (in several variants) and show that both of them are at least as hard as Partial MAP itself. We will also further elaborate on the exact complexity of decision versions of both problems which turns out to be non-trivial. We conclude in Section \ref{Conclusion}.

\section{Preliminaries}
\label{Preliminaries}

In this section we introduce our notational conventions. Specifically we will cover Bayesian networks, the complexity classes \ensuremath{\mathsf{PP}}\ and \ensuremath{\mathsf{NP^{\mathsf{PP}}}}, one-Turing reductions, and formal definitions of the approximation notions we will use in the paper. The reader is referred to textbooks like \cite{Darwiche09} (specifically complexity in Bayesian networks) and to \cite{Kwisthout15} (for a formal treatment of MAP approximations) for more background.

A Bayesian network $\ensuremath{\mathcal{B}} = (\mathbf{G}_{\mathcal{B}}, \ensuremath{\mathrm{Pr}})$ is a probabilistic graphical model that succinctly represents a joint probability distribution $\mprob{\mathbf{V}} = \prod_{i=1}^n \cprob{V_i}{\pi(V_i)}$ over a set of discrete random variables $\mathbf{V}$. $\ensuremath{\mathcal{B}}$ is defined by a directed acyclic graph $\mathbf{G}_{\mathcal{B}} = (\mathbf{V}, \mathbf{A})$, where $\mathbf{V}$ represents the stochastic variables and $\mathbf{A}$ models the conditional (in)dependences between them, and a set of parameter probabilities $\ensuremath{\mathrm{Pr}}$ in the form of conditional probability tables (CPTs). In our notation $\pi(V_i)$ denotes the set of parents of a node $V_i$ in $\mathbf{G}_{\mathcal{B}}$. We use upper case to indicate variables, lower case to indicate a specific value of a variable, and boldface to indicate sets of variables respectively joint value assignments to such a set.

One of the key computational problems in Bayesian networks is the problem to find the most probable explanation for a set of observations, i.e., a joint value assignment to a designated set of variables (the explanation set) that has maximum posterior probability given the observed variables (the joint value assignment to the evidence set) in the network. If the network includes variables that are neither observed nor to be explained (referred to as intermediate variables) this problem is typically referred to as \textsc{Partial MAP}. We use the following formal definition:

\fproblem{\textsc{Partial MAP}}
{A Bayesian network $\ensuremath{\mathcal{B}} = (\mathbf{G}_{\mathcal{B}}, \ensuremath{\mathrm{Pr}})$, where $\mathbf{V}(\mathbf{G}_{\mathcal{B}})$ is partitioned into a set of evidence nodes $\mathbf{E}$ with a joint value assignment $\mathbf{e}$, a set of intermediate nodes $\mathbf{I}$, and an explanation set $\mathbf{H}$}
{A joint value assignment $\mathbf{h}$ to $\mathbf{H}$ such that for all joint value assignments $\mathbf{h^{\prime}}$ to $\mathbf{H}$, $\cprob{\mathbf{h}}{\mathbf{e}} \geq \cprob{\mathbf{h^{\prime}}}{\mathbf{e}}$}

The following notation is taken from \cite{Kwisthout15}. For an arbitrary \textsc{Partial MAP}\ instance $\{\ensuremath{\mathcal{B}},\mathbf{H},\mathbf{E},\mathbf{I},\mathbf{e}\}$, let $\ensuremath{\mathit{cansol}_{\bn}}$ refer to the set of candidate solutions, with $\ensuremath{\mathit{optsol}_{\bn}}\in\ensuremath{\mathit{cansol}_{\bn}}$ denoting the {\em optimal} solution (or, in case of multiple solutions with the same posterior probability, one of the optimal solutions) to the \textsc{Partial MAP}\ instance; we will informally refer to this solution as the MAP explanation. When $\ensuremath{\mathit{cansol}_{\bn}}$ is ordered according to the probability of the candidate solutions (breaking ties between candidate solutions with the same probability arbitrarily), then $\ensuremath{\mathit{optsol}_{\bn}^{1 \dots m}}$ refers to the set of the first $m$ elements in $\ensuremath{\mathit{cansol}_{\bn}}$, viz. the $m$ most probable solutions to the \textsc{Partial MAP}\ instance.

We assume that the reader is familiar with standard notions in computational complexity theory, notably the classes \pt\ and \ensuremath{\mathsf{NP}}, \ensuremath{\mathsf{NP}}-hardness, and polynomial time (many-one) reductions. The class \ensuremath{\mathsf{PP}}\ is the class of decision problems that can be decided by a probabilistic Turing machine in polynomial time; that is, where \textit{Yes}-instances are accepted with probability strictly larger than $\sfrac{1}{2}$ and \textit{No}-instances are accepted with probability no more than $\sfrac{1}{2}$. A problem in \ensuremath{\mathsf{PP}}\ might be accepted with probability $\sfrac{1}{2} + \epsilon$ where $\epsilon$ may depend exponentially on the input size $n$. Hence, it may take exponential time to increase the probability of acceptance (by repetition of the computation and taking a majority decision) close to $1$. This is consistent with the sampling variant of the Chernoff bound: The number of samples $M$ needed to increase the probability of acceptance of \textit{Yes}-instances to $1-\delta$ is at least $\frac{\ln(\sfrac{1}{\sqrt{\delta}})}{\epsilon^2}$; when $\epsilon = \sfrac{1}{2^n}$ then $M$ is exponential in the input size. \ensuremath{\mathsf{PP}}\ hence is a powerful class; we know for example that $\ensuremath{\mathsf{NP}}\subseteq\ensuremath{\mathsf{PP}}$ and the inclusion is assumed to be strict. The canonical \ensuremath{\mathsf{PP}}-complete decision problem is \textsc{MajSAT}: given a Boolean formula $\phi$, does the majority of truth assignments to its variables satisfy $\phi$?

In computational complexity theory, so-called {\em oracles} are theoretical constructs that increase the power of a specific Turing machine. An oracle (e.g., an oracle for \ensuremath{\mathsf{PP}}-complete problems) can be seen as a `magic sub-routine' that answers class membership queries (e.g, in \ensuremath{\mathsf{PP}}) in a single time step. In this paper we are specifically interested in classes defined by non-deterministic Turing machines with access to a \ensuremath{\mathsf{PP}}-oracle. Such a machine is very powerful, and likewise problems that are complete for the corresponding complexity class \ensuremath{\mathsf{NP^{\mathsf{PP}}}}\ are highly intractable.

A decision variant of \textsc{Partial MAP}\ is known to be \ensuremath{\mathsf{NP^{\mathsf{PP}}}}-complete, even for binary variables, indegree at most $2$, and under the assumption that there exists at least one joint value assignment $\mathbf{h}$ such that $\mprob{\mathbf{h,e}} > 0$ \cite{Park04,Kwisthout11a}. In the intractability proofs in Section \ref{Main_results} we will assume, without loss of generality, that these constraints hold for all Bayesian networks $\ensuremath{\mathcal{B}}$ under consideration. As we use reductions from function problems, not decision problems, our reductions are formally polynomial-time one-Turing reductions. A function $f$ one-Turing reduces to $g$ (notation $f \leq^{\mathsf{FP}}_{1-T} g$) if there are functions $t_1$ and $t_2$ such that for all $x, f(x) = t_1(x, g(x, t_2(x)))$ \cite[p.5]{Toda94}.

We finish this section be repeating the following formal definition of rank-approximation of \textsc{Partial MAP}\ from \cite{Kwisthout15}; we will build on this definition in the next section.

\begin{definition}[rank-approximation of \textsc{Partial MAP}]
Let $\ensuremath{\mathit{optsol}_{\bn}^{1 \dots m}}\subseteq\ensuremath{\mathit{cansol}_{\bn}}$ be the set of the $m$ most probable solutions to a \textsc{Partial MAP}\ problem and let $\ensuremath{\mathit{optsol}_{\bn}}$ be the optimal solution. An explanation $\ensuremath{\mathit{approxsol}_{\bn}}\in\ensuremath{\mathit{cansol}_{\bn}}$ is defined to $m$-rank-approximate $\ensuremath{\mathit{optsol}_{\bn}}$ if $\ensuremath{\mathit{approxsol}_{\bn}}\in\ensuremath{\mathit{optsol}_{\bn}^{1 \dots m}}$.
\end{definition}

\section{Main results}
\label{Main_results}

Let $d_H$ be the Hamming distance between two joint value assignments. We define the following two problem variants to \textsc{Partial MAP}, where $m$ and $d$ are arbitrary constants:

\fproblem{$(m,d)$-\textsc{Similar Partial MAP}}
{As in \textsc{Partial MAP}}
{An explanation $\ensuremath{\mathit{approxsol}_{\bn}}\in\ensuremath{\mathit{cansol}_{\bn}}$ that $m$-rank-approximates $\ensuremath{\mathit{optsol}_{\bn}}$ and where $d_H(\ensuremath{\mathit{approxsol}_{\bn}}, \ensuremath{\mathit{optsol}_{\bn}}) \leq d$, or special symbol $\emptyset$ if no such explanation exists}

\fproblem{$(m,d)$-\textsc{Dissimilar Partial MAP}}
{As in \textsc{Partial MAP}}
{An explanation $\ensuremath{\mathit{approxsol}_{\bn}}\in\ensuremath{\mathit{cansol}_{\bn}}$ that $m$-rank-approximates $\ensuremath{\mathit{optsol}_{\bn}}$ and where $d_H(\ensuremath{\mathit{approxsol}_{\bn}}, \ensuremath{\mathit{optsol}_{\bn}}) \geq d$, or special symbol $\emptyset$ if no such explanation exists}

We will prove that both problems are \ensuremath{\mathsf{NP^{\mathsf{PP}}}}-hard by reduction from \textsc{Partial MAP}, even for binary variables with indegree at most $2$. We will start with the construction for $(m,d)$-\textsc{Similar Partial MAP}\ without the latter constraints and prove \ensuremath{\mathsf{NP^{\mathsf{PP}}}}-hardness (Figure \ref{construct}, panel a), and then adapt it to contain only binary variables and indegree at most $2$ (panel b). Then we show how this construction can be extended to prove \ensuremath{\mathsf{NP^{\mathsf{PP}}}}-hardness of $(m,d)$-\textsc{Dissimilar Partial MAP}\ (panel c).

\makefig{construct}
{Nodes added to \textsc{Partial MAP}\ network \ensuremath{\mathcal{B}}. a) basic construct; b) construct with binary nodes only; c) assuming at least $d$ different nodes.}
{10cm}{pdf}

\begin{theorem}
$(m,d)$-\textsc{Similar Partial MAP}\ is \ensuremath{\mathsf{NP^{\mathsf{PP}}}}-hard.
\label{spmap_thm}
\end{theorem}

\begin{proof}
We reduce from the \ensuremath{\mathsf{NP^{\mathsf{PP}}}}-hard \textsc{Partial MAP}\ problem. Let $\{\ensuremath{\mathcal{B}},\mathbf{H},\mathbf{E},\mathbf{I},\mathbf{e}\}$ be an instance to \textsc{Partial MAP}. From \ensuremath{\mathcal{B}}, we create an instance \ensuremath{\mathcal{B}'}\ to $(m,d)$-\textsc{Similar Partial MAP}\ as follows. To \ensuremath{\mathcal{B}}\ we add a singleton, unconnected node $M$ with $m$ uniformly distributed values $m_1 \ldots m_m$. Let $\mathbf{H'} = \mathbf{H} \cup \{M\}$ and observe that the MAP explanation $\mathbf{h}$ to $\mathbf{H}$ in the \textsc{Partial MAP}\ instance now translates into a set of $m$ best explanations $\ensuremath{\mathit{optsol}_{\bn}^{1 \dots m}} = \mathbf{h} \cup \{m_i\} (i = 1 \ldots m)$ with equal probability and by construction all these explanations differ only in the value assignment to $M$. Let $\ensuremath{\mathit{optsol}_{\bn}}$ be any arbitrary explanation in $\ensuremath{\mathit{optsol}_{\bn}^{1 \dots m}}$ and let $apr$ be a distinct explanation in $\ensuremath{\mathit{optsol}_{\bn}^{1 \dots m}}$. We have that $\ensuremath{\mathit{approxsol}_{\bn}}\in\ensuremath{\mathit{optsol}_{\bn}^{1 \dots m}}$ $m$-rank-approximates $\ensuremath{\mathit{optsol}_{\bn}}$ and that for any $d \geq 1$, $d_H(\ensuremath{\mathit{approxsol}_{\bn}}, \ensuremath{\mathit{optsol}_{\bn}}) = 1 \leq d$;  obviously the one-Turing reduction takes polynomial time and hence $(m,d)$-\textsc{Similar Partial MAP}\ is \ensuremath{\mathsf{NP^{\mathsf{PP}}}}-hard.
\end{proof}

Note that we can replace $M$ in this construction by $k = \lceil \log m \rceil$ uniformly distributed binary variables $M_i$ such that $\mathbf{H'} = \mathbf{H} \cup \{M_1 \ldots M_k \}$. A caveat here is that $\mathbf{h}$ may translate to more than $m$ explanations with equal probability and if $\ensuremath{\mathit{optsol}_{\bn}^{1 \dots m}}$ and $\ensuremath{\mathit{optsol}_{\bn}}$ are picked randomly from this set, we may end up\footnote{An example of such a situation would be when $m = 5$ and the solutions with binary encodings $000$, $011$, $101$, $110$, and $111$ would be in $\ensuremath{\mathit{optsol}_{\bn}^{1 \dots m}}$, with $\ensuremath{\mathit{optsol}_{\bn}} = 000$.} with a set that does not contain $\ensuremath{\mathit{approxsol}_{\bn}}$ such that $d_H(\ensuremath{\mathit{approxsol}_{\bn}}, \ensuremath{\mathit{optsol}_{\bn}}) = 1$. We therefore impose the constraint that $\ensuremath{\mathit{optsol}_{\bn}^{1 \dots m}}$ does not arbitrarily break ties, but that it contains the first $m$ explanations according to their lexicographical order, and because of this we can be sure that $\ensuremath{\mathit{optsol}_{\bn}^{1 \dots m}}$ contains at least one explanation $\ensuremath{\mathit{approxsol}_{\bn}}$ with $d_H(\ensuremath{\mathit{approxsol}_{\bn}}, \ensuremath{\mathit{optsol}_{\bn}}) = 1$.

\begin{corollary}
$(m,d)$-\textsc{Similar Partial MAP}\ is \ensuremath{\mathsf{NP^{\mathsf{PP}}}}-hard, even if all nodes are binary and have indegree at most 2.
\end{corollary}

We now extend the construction in the proof of theorem \ref{spmap_thm} to prove \ensuremath{\mathsf{NP^{\mathsf{PP}}}}-hardness of $(m,d)$-\textsc{Dissimilar Partial MAP}. To the above constructed network \ensuremath{\mathcal{B}'}\ we add $d-1$ nodes $D_1 \ldots D_{d-1}$ with $m$ uniformly distributed values. $D_1$ has $M$ as sole parent, whereas $D_i, i \geq 2$ has $D_{i-1}$ as sole parent. The probability of $D_i$ is defined to be deterministically dependent on its parent, i.e., $\cprob{D_1 = d_{1,j}}{M = m_j} = 1$ and $\cprob{D_i = d_{i,j}}{D_{i-1} = d_{i-1,j}} = 1$. We set $\mathbf{H'} = \mathbf{H} \cup \{M\} \cup \{D_1 \ldots D_{d-1}\}$. By construction, every explanation in $\ensuremath{\mathit{optsol}_{\bn}^{1 \dots m}}$ will differ in $d$ variables and thus any $\ensuremath{\mathit{approxsol}_{\bn}}\in\ensuremath{\mathit{optsol}_{\bn}^{1 \dots m}}$ $m$-rank-approximates $\ensuremath{\mathit{optsol}_{\bn}}$ and $d_H(\ensuremath{\mathit{approxsol}_{\bn}}, \ensuremath{\mathit{optsol}_{\bn}}) = d$. The construction with $k = \lceil \log m \rceil$ binary variables $M_i$ is similar, but now we add $(d-1)k$ nodes $D_{1,1} \ldots D_{d-1,k}$, which implies that every $\ensuremath{\mathit{approxsol}_{\bn}}\in\ensuremath{\mathit{optsol}_{\bn}^{1 \dots m}}$ differs in at least $d$ nodes.

\begin{corollary}
$(m,d)$-\textsc{Dissimilar Partial MAP}\ is \ensuremath{\mathsf{NP^{\mathsf{PP}}}}-hard, even if all nodes are binary and have indegree at most 2.
\end{corollary}

Note that $m$ and $d$ are constants and not part of the input. If we would make them part of the input, then unary notation would be necessary as the number of nodes added would otherwise be exponential in the binary representation of $m$ and $d$. However, since $m \leq 2^n$ (where $n$ denotes the number of variables in \ensuremath{\mathcal{B}}) the reduction would then not be polynomial in the \textsc{Partial MAP}\ instance any more. Hence, we require $m$ (and $d$) to be a constant.

\subsection{On membership in \ensuremath{\mathsf{NP^{\mathsf{PP}}}}}
\label{On_membership_in_nppp}

As explicated before, \textsc{Partial MAP}\ has an \ensuremath{\mathsf{NP^{\mathsf{PP}}}}-complete decision variant \cite{Park04}. Is it likely that appropriate decision variants of \textsc{Similar Partial MAP}\ and \textsc{Dissimilar Partial MAP}\ are also in \ensuremath{\mathsf{NP^{\mathsf{PP}}}}? We can check in polynomial time, using the MAP explanation \ensuremath{\mathit{optsol}_{\bn}}\ and the approximation $\ensuremath{\mathit{approxsol}_{\bn}}\in\ensuremath{\mathit{cansol}_{\bn}}$, that the Hamming distance $d_H(\ensuremath{\mathit{approxsol}_{\bn}}, \ensuremath{\mathit{optsol}_{\bn}})$ is as promised. However, a verification algorithm for $(m,d)$-\textsc{Similar Partial MAP}\ and $(m,d)$-\textsc{Dissimilar Partial MAP}\ should also verify that \ensuremath{\mathit{approxsol}_{\bn}}\ actually is within the $m$ best explanations, and finding the $m$-th best explanation is known to be \ensuremath{\mathsf{FP^{PP^{\mathsf{PP}}}}}-complete \cite{Kwisthout11a}. Hence, it is unlikely that $(m,d)$-\textsc{Similar Partial MAP}\ and $(m,d)$-\textsc{Dissimilar Partial MAP}\ have an \ensuremath{\mathsf{NP^{\mathsf{PP}}}}-complete decision variant.

\section{Conclusion}
\label{Conclusion}

In this short paper we elaborated on the computational complexity of finding MAP explanations that are almost-as-good as the most probable one (where `almost-as-good' was defined by rank) {\em and} that are sufficiently similar or dissimilar to the most probable explanation. We found that these problems are \ensuremath{\mathsf{NP^{\mathsf{PP}}}}-hard, that is, not easier than the \textsc{Partial MAP}\ problem itself. An attempt to extend this proof construct to include value-approximation did not succeed, as the proof constructs really require that all explanations in $\ensuremath{\mathit{approxsol}_{\bn}}\in\ensuremath{\mathit{optsol}_{\bn}^{1 \dots m}}$ have the same probability. In an extreme case it might be that we have just two explanations for a given \textsc{Partial MAP}\ problem, where $\cprob{h}{\mathbf{e}} = \sfrac{1}{2} + \epsilon$ and $\cprob{\neg h}{\mathbf{e}} = \sfrac{1}{2} - \epsilon$, and any non-uniform probability assignment to the additional variables in the explanation set will destroy the proof.

In future work we'd like to extend our results to value-approximation as well as rank-approximation. It might be relevant to investigate the complexity of these problems in constrained structures, such as polytrees, where the \textsc{Partial MAP}\ problem nonetheless remains \ensuremath{\mathsf{NP}}-hard.

\bibliographystyle{splncs04}


\section{Introduction}

One of the major unsolved issues in astrophysics is how galaxy buildup proceeded in the early Universe. Evolutionary models predict that different mechanisms can play a fundamental role for galaxy growth, such as galaxy mergers \citep[e.g.][]{hopkins2006, somerville2008} or cold gas accretion within gas-rich proto-disks \citep[e.g.][]{dekel2009}. Whatever the physical mechanisms dominating galaxy evolution might be, they must reproduce the observed distribution of baryons at high redshift, and connect it to the subsequent evolution of galaxies within dark matter halos. In this context, it becomes crucial to have accurate measurements of galaxy number densities and stellar-mass evolution over most of the cosmic history ($\mathrm{z < 3 - 4}$). So far, important discrepancies still exist between predictions from theoretical models and observations of the number density of the most massive galaxies ($\mathrm{M_{\star}>3×10^{11} M_{\odot}}$) at high redshifts \citep[eg.,][]{fontanot2009, marchesini2009,marchesini2010, ilbert2013, hirschmann2016}. In fact, theoretical models still under-predict the observed number density of these galaxies, unless the model-predicted stellar mass functions (SMFs) are convolved with large errors in  $\mathrm{M_{\star}}$ \citep[e.g., $\sim \, 0.25-0.3$ dex at z=2-3,][]{marchesini2009,marchesini2010,ilbert2013, muzzin2013b}, arguably much larger than the typical random errors on $\mathrm{M_{\star}}$ for massive (hence bright) galaxies.
This disagreement at the very high-mass end has been recently confirmed by recent surveys, e.g. the NEWFIRM (NOAO Extremely Wide-Field Infrared Imager ) Medium-Band Survey (NMBS) \citep{whitaker2011} and the UltraVISTA survey \citep{muzzin2013a, hirschmann2016}.
NMBS is a moderately deep survey covering 0.45 $\mathrm{deg^2}$ over two $\mathrm{30\arcmin \times 30 \arcmin}$ fields within the Cosmic Evolution Survey \citep[COSMOS;][]{scoville2007} and the All-Wavelength Extended Groth Strip International Survey \citep[AEGIS;][]{davis2007}, chosen to take advantage of the wealth of publicly available ancillary data over a broad wavelength range, from the X-ray to the radio.
 Thanks to the medium-band near-infrared (NIR) filters,  NMBS provides precise $\mathrm{z_{phot}}$ ($\mathrm{\Delta z/(1 + z)\approx 0.02}$) and well-sampled spectral energy distributions (SEDs) for thousands of K-selected galaxies at $z>1.5$. \\
One of the most remarkable results from NMBS is the clear detection of the galaxy bi-modality in the galaxy population (quiescent vs. star-forming) out to $\mathrm{z \sim 2.5}$ \citep{brammer2009}.
The NMBS showed that nearly all galaxies with $\mathrm{log(M/M_{star})\, > \, 10.5 }$ have red rest-frame U−V colors up to z $\mathrm{=}$ 2.  Despite similar U−V colors, the NMBS data showed that these galaxies have strikingly different properties, with approximately half being red due to a lack of recent star formation (quiescent), the other half being red due to vigorous dusty star formation. This result showed the limitations of using a single color to identify quiescent galaxies at $\mathrm{z\,>\,1}$. Instead, star-forming and quiescent galaxies can be robustly separated out to $\mathrm{z \sim\, 3}$ with a two-color criterion using the rest-frame U$-$V versus V$-$J diagram (UVJ diagram, hereafter).
However, for high redshift galaxies, $\mathrm{3.6\, and\,  4.5\, \mu m}$ imaging with IRAC \citep[][]{fazio2004} aboard the Spitzer Space Telescope \citep{werner2004} are absolutely critical for the UVJ diagnostic, as the rest-frame V $-$ J color cannot be directly measured without IRAC for galaxies at $\mathrm{0.8\,< \, z\, <\, 3}$.
It is only thanks to IRAC in combination with the very precise photometric redshift ($\mathrm{z_{phot}}$) from NMBS that the existence of the bi-modality and the existence of quiescent galaxies out to z$\sim$2.5 (and beyond) have been shown definitively \citep{whitaker2011}. The lack of IRAC photometry prevents robust classification into quiescent and star-forming galaxies at $z>0.8$.\\
While our knowledge of massive galaxies has grown significantly in the past years due to NMBS and other many NIR surveys, little is known about the evolution of the most massive galaxies  (UMGs) with $\mathrm{M_{\star} >3\times10^{11}M_{\odot}}$.  This is simply because their space density is extremely low ($\sim$30× lower than galaxies with $\mathrm{M_{\star} =10^{11}M_{\odot}}$). Even the largest NIR surveys have lacked the volume to detect them in significant numbers. Despite their rareness, these galaxies arguably represent the biggest problem in our understanding of how galaxies form. Abundance matching models suggest that  ultra-massive galaxies, those with $\mathrm{M_{\star} =3\times 10^{11}M_{\odot}}$  should reside in dark matter haloes of a few $\mathrm{\times 10^{14}M_{\odot}}$ at all redshifts, i.e., they are the progenitors of local Brightest Cluster Galaxies. While observationally the most massive galaxies appear to be the first to become quiescent, it is completely unknown why this is the case and what drives the quenching of the star formation. \\
In 2010-2012, NMBS-II (PI: Van Dokkum),  a 44-night NOAO survey program with NEWFIRM was conducted with the primary goal of surveying a large enough volume to identify a statistically robust sample of UMGs out to $z=3.5$. NMBS-II is shallower but $11\times$ wider than NMBS, and it is unique for its combination of NIR medium-band imaging, well-sampled SEDs, and large area. NMBS-II covers a total of $\sim$ 5.4 $\mathrm{deg^2}$ over 7 fields, namely CFHTLS-D1, CFHTLS-D2, (COSMOS), CFHTLS-D3, CFHTLS-D4, and 3 MUSYC (Multiwavelength Survey by Yale-Chile) fields, namely MUSYC1030, MUSYC1255 and MUSYC-ECDF.
Observations of the NMBS-II are now completed, and multi-wavelength K-selected catalogs have been constructed (including deep optical ugriz data from the CFHTLS and U-V data from Galex). The NMBS-II J1J2J3H1H2K data, combined with deep optical data from the CFHTLS and IRAC imaging over the full survey area, promise to provide a mass-complete sample of $\sim$ 300 UMGs ($\mathrm{log(M_{\star}/M_{\odot}) > 11.4}$, of which $\sim$50 with $\mathrm{log(M_{\star}/M_{\odot}) > 11.6}$) with high-confidence $\mathrm{z_{phot}}$ at $1.5<z<3$. This sample, spanning the full range of star-formation activity, will provide the first comprehensive view of the population of UMGs at $z>1.5$. This will be the largest sample of UMGs with accurate $\mathrm{z_{phot}}$ ever constructed.
NMBS-II is truly unique as no other current survey has its combination of surveyed area and NIR medium-band coverage, preventing the delivery of sufficiently precise $\mathrm{z_{phot}}$ over a wide enough area to assemble similarly large samples of UMGs at $z>1.5$.\\
However, Spitzer-IRAC coverage from previous programs of NMBS-II fields is available only for $\sim$60\% of the surveyed  area, hence missing for $\sim 2.2 \mathrm{deg^2}$.  IRAC data, besides being critical to classify the population of galaxies at z$>$0.8 into quiescent and star-forming galaxies using the UVJ diagram, are crucial for more accurate estimates of $z_{phot}$ and $M_{\star}$, as well as for robust determinations of SFRs, as they significantly help breaking the age-dust degeneracy in the modeling of the SEDs.   \\
With this paper, we present and publicly release the scientific mosaics, coverage maps, and 3$\sigma$ depth maps  at 3.6$\mu$m and 4.5$\mu$m for five of the NMBS-II fields. The mosaics combine both archival data and new data observed in Cycle 10.
The paper is organized as follows. Section \ref{s:obs} describes all the programs that have been used to create single, contiguous deep images in the 3.6 and 4.5$\mu m$ bands. In Section \ref{s:reduction}, we describe how the mosaics and coverage maps were produced, and Section~\ref{s:prop} presents the characterization of the mosaics, including image quality and depth.
All coordinates refer to the J2000 system. All magnitudes are in the AB system.

\section{ Observations}\label{s:obs}
In this section,  we present the NMBS-II fields and the individual Spitzer-IRAC programs that have been combined to produce deep IRAC mosaics over the NMBS-II, completing coverage of the full survey. Although the NMBS-II survey is comprised of seven fields, we produce IRAC mosaics only for five fields, because comprehensive IRAC mosaics in COSMOS and ECDFS are already publicly available from the Spitzer Science Center as part of previous IRAC surveys.

\subsection{NMBS-II Fields}\label{ss:nmbs-ii fields}
In this sub-section, we present the seven NMBS-II fields, including the ancillary data publicly available.

\begin{itemize}
\item CFHTLS-D1 (XMM-LSS) is a $\mathrm{1 \, deg^2}$ field centered on RA= 02:26:00.00 (hh:mm:ss) and DEC= -04:30:00.00 (dd:mm:ss). This field has ugriz imaging from the Canada-France-Hawaii Telescope Legacy Survey (CFHTLS).
 This field has NIR  images in YJHK taken  as part of the Visible and Infrared Survey Telescope for
Astronomy (VISTA) Deep Extragalactic Observations (VIDEO) survey \citep{jarvis2013}. Other NIR data are available through the WIRCam Infrared Deep Survey \citep[WIRDS;][]{bielby2012}. Imaging of the CFHTLS-D1 field at 250, 350, and 500 $\mu$m was
provided by the Spectral and Photometric Imaging REceiver \citep[SPIRE;][]{griffin2010} aboard the Herschel Space Observatory \citep{pilbratt2010}, taken as part of the Herschel Multi-tiered Extragalactic Survey \citep[HerMES;][]{oliver2012} X-ray observations of the CFHTLS-D1 field were taken from the X-ray Multi-mirror Mission space telescope \citep[XMM-Newton;][]{jansen2001} as part of the XMM Medium Deep Survey \citep[XMDS;][]{chiappetti2005}. Spectroscopic data for this field have been acquired with the VIsible Multi-Object Spectrograph (VIMOS) as part of the VIMOS-VLT deep survey (VVDS).
As a part of NMBS-II, this field was imaged in medium-bandwidth filters ($J_1J_2J_3H_1H_2$) over the wavelength range $\mathrm{1–1.8\mu m}$ and in the K-band.
CFHTLS-D1 has a  100\% data coverage from the Spitzer Wide-area InfraRed Extragalactic Survey \citep[SWIRE,][]{lonsdale2003}. SWIRE is a wide-area, high galactic latitude imaging survey which covered  roughly 50 $\mathrm{deg^2}$ across a number of fields with the MIPS far-infrared camera and the with the IRAC mid-infrared camera.

\item CFHTLS-D2  is a $\mathrm{2 \, deg^2}$ field centered on  RA= 10:00:28.60 and DEC= 02:12:21.00. This field is part of the Cosmic Evolution Survey survey. \citep[COSMOS;][]{scoville2007}. Within COSMOS, $\mathrm{1\, deg^2}$ has is also been imaged by CFHT as part of the Canada-France-Hawaii Telescope Legacy Survey Deep (CFHTLS-Deep). This field has X-ray imaging \citep[Chandra and XMM,][respectively]{elvis2009,finoguenov2007} from COSMOS, HST imaging from the original COSMOS Hubble Space Telescope Treasury project, and from the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey \citep[CANDELS;][]{koekemoer2011}, deep optical imaging taken with Hyper Suprime-Cam \citep[HSC;][]{tanaka2017}, and deep NIR imaging as part of UltraVISTA \citep{mckracken2013} survey.  As a part of NMBS-II, this field was imaged in medium-bandwidth filters ($\mathrm{J_1 J_2 J_3 H_1 H_2}$),  and in the K-band.
Narrow-band NIR images of the COSMOS-HST field we also taken with Suprime-Cam on the Subaru Telescope within the Subaru COSMOS 20 project \citep{taniguchi2015}. COSMOS has extensive spectroscopic data with over 97,000 spectra targeting over 68,000 unique objects.  Major surveys include zCOSMOS \citep{lilly2007}, and the VIMOS Ultra Deep Survey \citep[VUDS;][]{lefevre2015}. HST spectroscopy is also available from the 3D-HST program \citep{vandokkum2011}. CFHTLS-D2 has 100\% IRAC data coverage coming from multiple surveys, e.g. S-COSMOS \citep{sanders2007}, the Spitzer Extended Deep Survey  \citep[SEDS;]{ashby2013}, S-CANDELS, \citep{ashby2015}, Star Formation at 4 $<$ z $<$ 6 from the Spitzer Large Area Survey with Hyper-Suprime-Cam \citep[SPLASH;][]{steinhardt2014}, and Spitzer Matching survey of the UltraVISTA ultra-deep Stripes \citep[SMUVS;][]{ashby2018}.

\item CFHTLS-D3  is a $\mathrm{1\, deg^2}$ field centered on RA= 14:19:27.00 and DEC= +52:40:56.00. This field was imaged by CFHT both as part of the AEGIS survey  and CFHTLSD. X-ray images of this field were taken with Chandra as part of the  AEGIS-X \citep{laird2009}.  Deep HST images of the EGS were obtained with the Advanced Camera for Surveys (ACS) as part of GO program 10134 (PI: M. Davis). As a part of NMBS-II, this field was imaged in medium-bandwidth filters ($\mathrm{J_1 J_2 J_3 H_1 H_2}$),  and in the K-band. CFHTLS-D3 previously had 20\% IRAC data coverage from  S-CANDELS, \citep{ashby2015}.

\item CFHTLS-D4  is a $\mathrm{1 \, deg^2}$ field centered on RA= 22:15:31.00  and DEC= -17:43:56.00. This field was imaged by CFHT as part of CFHTLS. As a part of NMBS-II, this field was imaged in medium-bandwidth filters ($\mathrm{J_1 J_2 J_3 H_1 H_2}$),  and in the K-band.
CFHTLS-D4 had no IRAC coverage previous to this paper.

\item MUSYC Extended CDFS (MUSYC-ECDFS) is a 0.5 $\mathrm{deg^2}$ field centered on RA= 03:32:29.0 and DEC= -27:48:47.00. This is one of the  fields of the Great Observatories Origins Deep Survey \citep[GOODS;][]{dickinson2001}. Within GOODS, the ACS camera on
HST acquired imaging of the ECDFS field. Ground-based optical imaging of the ECDFS with the Wide Field Imager (WFI) on the ESO/MPG 2.2m telescope was obtained as part of the COMBO-17 \citep{wolf2004}. NIR images of this field were also taken with the VLT/ISAAC.
This field has also NIR  images in YJHK taken  as part of VIDEO survey \citep{jarvis2013}. As a part of NMBS-II this field was imaged in medium-bandwidth filters ($\mathrm{J_1 J_2 J_3 H_1 H_2}$),  and in the K-band.
MUSYC-ECDFS has 100\% data coverage from the Spitzer IRAC/MUSYC Public Legacy Survey in the Extended Chandra Deep Field South \citep[SIMPLE,][]{damen2011}. SIMPLE consists of deep IRAC observations (several hours per pointing) covering the 0.5 x 0.5 $\mathrm{deg^2}$ area surrounding the GOODS CDFS field.
\item MUSYC1030 is a 0.5 $\mathrm{deg^2}$ field centered around RA= 10:30:27.10 and DEC= 05:24:55.02. The center of this field is the position of the SDSS quasar SDSS J1030+0524S \citep{becker2001}. This field was imaged by different telescopes within several surveys. Optical images of this field were taken within the MUSYC survey \citep{gawiser2006a}. Complimentary K-band images were taken using the Infrared Side Port Imager (ISPI) on the 4m Blanco Telescope at the Cerro Tololo Inter-American Observatory (CTIO) in the MUSYC wide NIR survey. As a part of NMBS-II, this field was imaged in medium-bandwidth filters ($\mathrm{J_1 J_2 J_3 H_1 H_2}$),  and in the K-band. MUSYC1030 has 10\%  IRAC coverage.

\item MUSYC1255 is a 0.5 $\mathrm{deg^2}$ centered at RA= 12:55:40.0 and DEC= 01:07:00.00 and also has optical and NIR  imaging as part of the
MUSYC and MUSYC-NIR surveys. As a part of NMBS-II, this field was imaged in medium-bandwidth filters ($\mathrm{J_1 J_2 J_3 H_1 H_2}$),  and in the K-band. MUSYC1255 has 10\%  IRAC coverage from the Spitzer Space Telescope Cycle-3 program GO-30873 (PI: Labbe ).
\end{itemize}

\subsection{New IRAC imaging in NMBS-II}\label{ss: new_irac}
To fully exploit the NMBS-II,  in Cycle 10, we were awarded 22 hours on Spitzer to complete the IRAC coverage of the NMBS-II fields in $\mathrm{3.6\, \mu m}$ (3.6$\mu$m) and $\mathrm{4.5\, \mu m}$ (4.5$\mu$m) (ID:10084; PI: D. Marchesini). The fields imaged in this program were: CFHTLS-D3, CFHTLS-D4,MUSYC1030, andMUSYC1255. The total area mapped in $\mathrm{3.6\, \mu m}$ and $\mathrm{4.5\, \mu m}$ in Cycle 10  is $\sim$ 2.2 $\mathrm{deg^2}$, for a total of 342 individual pointings.\\
In this paper, we present scientific mosaics and coverage maps for CFHTLS-D1,  CFHTLS-D3, CFHTLS-D4,MUSYC1030, andMUSYC1255.

\section{Data reduction}\label{s:reduction}

We downloaded all archival data from the Spitzer Heritage Archive (SHA) and combined them with our data sets. We reduced all data in a consistent manner, coadding them into one mosaic per field. We followed the procedure described in \citep{labbe2015}.
The individual exposures are organized into self-contained segments known as Astronomical Observing Requests (AORs). Over the five fields, there are 755 AORs for each filter.  Each AOR consists of basic calibrated data (BCD). BCDs are exposure-level data after having passed through the pipelines. Instrumental signatures have been mostly removed, and the BCDs have been absolutely calibrated into physical units (i.e., $MJy/sr =  10^{-17} erg s^{-1} cm^{-2} Hz^{-1} sr^{-1}$).  The BCDs were improved by applying a series of procedures in order to remove cosmic rays and remaining instrumental artifacts including column pulldown,  muxbleed, maxstripe, and the ``first-frame effect'' \citep{hora2004}.  Single background-subtracted frames were then combined into a mosaic image large enough to hold all input frames (using WCS astrometry  to align the images). In this step,  we masked bad pixels and applied a distortion correction in WCS. This image, in combination with the reference image, was then used to refine the pointing of the individual mosaics. The individual pointing-refined frames were registered to and projected on the reference image. \\
The reference images for the CFHTLSD-1, CFHTLSD-3, CFHTLSD-4 fields are the $z-$band stack from CFHTLSDeep taken from the CFHTLS  T0007 release. These stacks are the combinations of the 25\% best seeing images (seeing $< 0.95 \arcsec$) of each fields in the  $z-$band.  Each stack comprises all images located inside a radius of $3 \arcmin$ with respect to a tile center position,  weighted accordingly, and combined using SWarp, with a Lanczos3 interpolation kernel. The total exposure times for each of the reference images is 14.20hr, 12.51hr and 12.00 hr, while the seeing is 0.559$\arcsec$, 0.537 $\arcsec$, and 0.551$\arcsec$, for CFHTLS-D1, CHFTLS-D3, and CFHTLS-D4 respectively. The pixel scale is 0.186 $\arcsec / pixel$. These images are available on the Terapix database (\url{http://terapix.iap.fr/cplt/T0007/table_syn_T0007.html}).\\
The final IRAC mosaics and coverage maps are then re-sampled to a higher pixel scale value which is listed in Table~\ref{tab-appiracdata1}.
The reference images for MUSYC1030 and MUSYC1255 are taken with  the Infrared Sideport Imager (ISPI)  camera mounted at the prime focus of the Blanco 4-m telescope at CTIO in the K-band filter. The original detector pixel scale is 0.305$\arcsec$/pixel. However both K-band images have been previously processed so that the pixel scale is changed slightly to 0.267 and 0.269 $\arcsec$/pixel, respectively for MUSYC1030 and MUSYC1255.
Table~\ref{tab-appiracdata1} summarizes the characteristics of the IRAC data over the five observed fields, such as the total exposure time, the number of AORs that contributed to the mosaic, and the pixel scale at $\mathrm{3.5 \, \mu m}$ and $\mathrm{4.5 \, \mu m}$. For comparison, we also add, where available, the characteristics of the other two fields.\\

\begin{deluxetable*}{ccccccc}
\tablewidth{0pt}
\tablecaption{Characteristics of the IRAC observations
		\label{tab-appiracdata1}}
\tablehead{  \colhead{Filter} & \colhead{PIDs\tablenotemark{a}}  & \colhead{Max Exposure Time\tablenotemark{b}} & \colhead{Number of AORs}  &  \colhead{Pixel scale} & \colhead{ZP} \\
	    	  ($\mu$m) & & (h) & & (arcsec/pixel) & (AB) }
 \startdata
 CFHTLS-D1: & & &  & &\\
3.6 & 181 & \textnormal{8.18} & 262  & 0.558 & 20.04\\
4.5 & 181 & \textnormal{7.31} & 262  & 0.558 & 20.04\\
 CFHTLS-D2\tablenotemark{c}: & & &  & &\\
3.6 & 20070, 61043,80057, 90042, 10159, 11016 &  &   & 0.6 & \\
4.5 & 20070, 61043,80057, 90042, 10159, 11016 &  &   & 0.6 & \\
\hline
CFHTLS-D3: & &  & & &\\
3.6 & 8, 41023, 61042, 80216, 10084 & \textnormal{29.3} & 456 & 0.558 & 20.04\\
4.5 & 8, 41023, 61042, 80216, 10084 & \textnormal{26.8} & 456  & 0.558 & 20.04 \\
\hline
CFHTLS-D4: & &  & & &\\
3.6 & 10084  & \textnormal{1.32} & 13  & 0.558 & 20.04 \\
4.5 & 10084  & \textnormal{0.89}  & 13  & 0.558 & 20.04\\
\hline
MUSYC-ECDF\tablenotemark{d}: & &  & & &\\
3.6 & 20708 & \textnormal{2.5} &   & 0.6 & 22.42\\
4.5 & 20708 & \textnormal{2.5} &   & 0.6 & 22.19\\
\hline
MUSYC1030: & &  & & &\\
3.6 & 3198, 10084 & \textnormal{1.70} & 11  & 0.267 & 20.04\\
4.5 & 3198, 10084 & \textnormal{1.76} & 11  & 0.267 & 20.04\\
\hline
MUSYC1255:&  &  & & &\\
3.6 & 30873, 10084 & \textnormal{1.40} & 13  & 0.269 & 20.04 \\
4.5 & 30873, 10084 & \textnormal{1.46} & 13  & 0.269 & 20.04\\
\hline
 \enddata
  \tablenotetext{a} {Spitzer Program Identification Numbers}
 \tablenotetext{b}{Maximum exposure time per position on sky per channel.}
 \tablenotetext{c}{The characteristics of this field are taken from \citealt{ashby2018}}
  \tablenotetext{d}{The characteristics of this field are taken from \citealt{damen2011}}
\end{deluxetable*}
\noindent In Figures \ref{f:cfhtlsd1}, \ref{f:cfhtlsd3}, \ref{f:cfhtlsd4}, \ref{f:musyc1030}, and \ref{f:musyc1255}, we show the reference image, the two IRAC mosaics at 3.6 $\mu  m$ and 4.5 $\mu  m$ and the corresponding coverage maps for each of the five fields.

\begin{figure*}[ht]
	\centering
	\includegraphics[scale=0.5]{./figures_journal/cfhtlsd1.pdf}
	\caption{ CFHTLS-D1: reference band mosaic, 3.6$\mu m$ and 4.5$\mu m$ mosaics and coverage maps. The coverage maps show where the individual pointings of IRAC overlap when they are combined to form the mosaic. Pixels with darker colors have more coverage.}
	\label{f:cfhtlsd1}
\end{figure*}

\begin{figure*}[ht]
	\centering
	\includegraphics[scale=0.5]{./figures_journal/cfhtlsd3.pdf}
	\caption{Same as in Figure~\ref{f:cfhtlsd1}, but for the CFHTLS-D3 field. For this field, we use a logarithmic scale for the coverage maps to show simultaneously the regions with the highest and lowest coverage.}
	\label{f:cfhtlsd3}
\end{figure*}

\begin{figure*}[ht]
	\centering
	\includegraphics[scale=0.5]{./figures_journal/cfhtlsd4.pdf}
	\caption{ Same as in Figure~\ref{f:cfhtlsd1}, but for the CFHTLS-D4 field. }
	\label{f:cfhtlsd4}
\end{figure*}

\begin{figure*}[ht]
	\centering
	\includegraphics[scale=0.5]{./figures_journal/musyc1030.pdf}
	\caption{ Same as in Figure~\ref{f:cfhtlsd1}, but for the MUSYC1030 field.}
	\label{f:musyc1030}
\end{figure*}

\begin{figure*}[ht]
	\centering
	\includegraphics[scale=0.5]{./figures_journal/musyc1255.pdf}
	\caption{ Same as in Figure~\ref{f:cfhtlsd1}, but for the MUSYC1255 field. }
	\label{f:musyc1255}
\end{figure*}

\section{Mosaic quality}\label{s:prop}
In this section, we describe the quality of the final mosaics  in both channels for  each of the five NMBS-II fields.

\textnormal{To characterize the FWHM and the PSF of these mosaics we proceed in this way. First,  we define a sample  of $\sim$ 50 bright, isolated and unsaturated stars in each channel and each field. Then, for each star in the sample we measured a curve of growth,  i.e. the distribution of the light at different radii, normalized by the flux of the star in aperture of sufficiently large diameter ($\mathrm{6 \arcsec}$). For each field,  we consider the median of this growth curves in each channel (black curves in Fig.~\ref{f:psf_visual_1} and Fig.~\ref{f:psf_visual_2}). Outliers, such as saturated stars, were then determined based on the shape of their light profile compared with the median curve of growth, and rejected from the sample.  We use the median growth curve in each field to derive the median half-light and 75\%-light radii, i.e.,  the radii that contain 50\% and 75\% of the light,  respectively. These values are reported in Table  ~\ref{t:psf}.  The same sample of stars is also used to derive the FWHMs of the IRAC images. The derived FWHM are listed in Table ~\ref{t:psf}. We stress that the IRAC PSF is complicated and that spatial variations across the field of view are to be expected given we mosaiced IRAC data taken from several independent programs. However, the state-of-the-art approach to extract IRAC photometry from wide and deep surveys consists in the adoption of software that accounts for spatially varying PSF across the field of view, such as MOPHONGO \citep{labbe2006}, TFIT \citep{laidler2007}, and T-PHOT \citep{merlin2015}. }
\begin{figure*}[h]
	\centering
	\includegraphics[scale=0.45]{./figures_journal/psf_visual_ch1.pdf}
	\caption{ First five panels: Growth curves for a sample of $\sim$ 50 bright unsaturated stars for all five fields at 3.6$\mathrm{\mu m}$ normalized to the flux of the star in an aperture of $\mathrm{6 \arcsec}$ diameter. The black curve in each panel represents the median of the distribution. The bottom-right panel shows the median growth curve of each field at  3.6$\mathrm{\mu m}$. }
	\label{f:psf_visual_1}
\end{figure*}

\begin{figure*}[h]
	\centering
	\includegraphics[scale=0.45]{./figures_journal/psf_visual_ch2.pdf}
	\caption{ Same as in Fig. ~\ref{f:psf_visual_1}, but at 4.5$\mathrm{\mu m}$. }
	\label{f:psf_visual_2}
\end{figure*}
In Table ~\ref{t:ap}, we list the  3$\sigma$ magnitude depth for all mosaics.  For comparison, we also add, where available, the FWHMs and  3$\sigma$ magnitude depths of the other two fields taken from the literature papers mentioned in Sect.~\ref{s:obs}. \\
\noindent We use the ``empty aperture'' method to empirically determine the noise properties of our IRAC mosaics, following the same approach as for the deep NIR MUSYC data \citep{quadri2007}.
Briefly, we randomly place a large number of apertures on the noise-normalized images (obtained by multiplying the images by the square root of the coverage maps), reject all apertures falling on sources, and measure the flux in the remaining ones. We chose an aperture size of $3\arcsec$ and measure the flux in $\sim$ 8000 apertures for CFHTLSD fields and  $\sim$ 4000 in the MUSYC fields.  A histogram of the measured fluxes is constructed for each field  in both $\mathrm{3.6\, \mu m}$ and $\mathrm{4.5\, \mu m}$, and shown in Figures~\ref{f:empty_cfhtlsd1}, ~\ref{f:empty_cfhtlsd3}, ~\ref{f:empty_cfhtlsd4}, ~\ref{f:empty_musyc1030} and ~\ref{f:empty_musyc1255} (bottom panels).
These histograms are all well fitted by a Gaussian.  The top panels of Figures ~\ref{f:empty_cfhtlsd1}, ~\ref{f:empty_cfhtlsd3}, ~\ref{f:empty_cfhtlsd4}, ~\ref{f:empty_musyc1030} and ~\ref{f:empty_musyc1255} show  the maps of  $\mathrm{3\sigma}$ magnitude depths derived using the empty aperture method according to the Equation:

\begin{equation}
\mathrm{depth(3\sigma)[AB] = -2.5 \log{(\frac{3\sigma_{best-fit}}{\sqrt{\mathrm{COVERAGE\, MAP}}})} + ZP,}
\end{equation}

where ZP is listed in Table~\ref{tab-appiracdata1} for each field.
Table ~\ref{t:ap} lists the 15th, 50th, and 75th percentiles of the 3$\sigma$ magnitude depths in both $\mathrm{3.6\, \mu m}$ and $\mathrm{4.5\, \mu m}$ in an aperture of D=3$\arcsec$ as derived with the empty aperture method. The median value of the 3$\sigma$ depths for $\mathrm{3.6\, \mu m}$ varies from 22.83 mag (in CFHTLSD4 ) to 24.52 (in CFHTLSD1), while for $\mathrm{4.5\, \mu m}$ it varies from 22.85 mag (in CFHTLSD4) to 24.44 (in CFHTLSD1 ).
As expected, regions with higher exposure times (see Figures~\ref{f:cfhtlsd1}, ~\ref{f:cfhtlsd3}, ~\ref{f:cfhtlsd4},~\ref{f:musyc1030}, and \ref{f:musyc1255})  also have fainter $3\sigma$ magnitude depths. As shown in Figure~\ref{f:empty_cfhtlsd3}, using both the images observed in this program and the data collected by the AEGIS survey, we present results in an extremely deep mosaic in some of the regions of CFHTLSD3.  \\

\begin{deluxetable*}{cccc}
	\tablewidth{0pt}
	\tablecaption{Characteristics of the IRAC observations
		\label{t:psf}}
	\tablehead{\colhead{Filter} & \colhead{FWHM} & \colhead{50\% light radius}  & \colhead{75\% light radius}  \\
		($\mu$m)      & (\arcsec)      & ( \arcsec)   & ( \arcsec) }
	\startdata
	CFHTLS-D1:& & &   \\
	3.6 &  2.08 & 2.20 & 3.28  \\
	4.5 &  1.95 & 2.20 & 3.39  \\
	\hline
	CFHTLS-D2\tablenotemark{a}:& & &  \\
	3.6 &  $\mathrm{\sim}$ 2.00 &  &   \\
	4.5 &  $\mathrm{\sim}$ 2.00 &  &    \\
	\hline
	CFHTLS-D3:& & & \\
	3.6 & 2.00 & 2.15 & 3.22\\
	4.5 & 1.86 & 2.16 & 3.38  \\
	\hline
	CFHTLS-D4:& & & \\
	3.6 & 1.86 & 2.08 & 3.17  \\
	4.5 & 1.75 & 2.14 & 3.35 \\
	\hline
	MUSYC-ECDF\tablenotemark{b}:& & &\\
	3.6 & 1.97 & &\\
	4.5 & 1.93 & & \\
	\hline
	MUSYC1030:& & &\\
	3.6 & 2.17 & 2.11 & 3.18  \\
	4.5 & 2.01 & 2.16 & 3.35  \\
	\hline
	MUSYC1255:& & &\\
	3.6 & 2.02 & 2.09 & 3.17 \\
	4.5 & 1.98 & 2.11 & 3.34  \\
	\hline
	\enddata
	\tablenotetext{a}{The FWHM of this field  are taken from \citealt{ashby2018}.}
	\tablenotetext{b}{The FWHM of this field are taken from \citealt{damen2011}. }
\end{deluxetable*}

\begin{deluxetable*}{cccc}
	\tablewidth{0pt}
	\tablecaption{Depths of the IRAC observations
		\label{t:ap}}
	\tablehead{\colhead{Filter} & \colhead{15th percentile of 3$\sigma$ depth\tablenotemark{a}}  & \colhead{Median of 3$\sigma$ depth \tablenotemark{a}} & \colhead{75th percentile of 3$\sigma$ depth\tablenotemark{a}} \\
		($\mu$m)          & ( AB)   & ( AB)  & ( AB)}
	\startdata
	CFHTLS-D1:& & &  \\
	3.6 &   24.33 & 24.52 & 24.66  \\
	4.5 &   24.24 & 24.44 &  24.58  \\
	\hline
	CFHTLS-D2\tablenotemark{b}:&  & &  \\
	3.6 &    & 25.31 &   \\
	4.5 &    & 25.31  &   \\
	\hline
	CFHTLS-D3:& & &\\
	3.6 & 22.79 & 23.76 &  25.54\\
	4.5 &  22.88 & 23.45 & 25.29  \\
	\hline
	CFHTLS-D4:&  & & \\
	3.6 &  22.53 & 22.83 & 23.00 \\
	4.5 &  22.51 & 22.85 &  23.03 \\
	\hline
	MUSYC-ECDF\tablenotemark{c}:& & &\\
	3.6 & 24.21  & 24.41 & 24.55 \\
	4.5 &  24.05 & 24.24 &  24.37 \\
	\hline
	MUSYC1030:&  & &\\
	3.6 &  22.63 & 22.92 & 23.20 \\
	4.5 &  22.63 & 22.97 &  23.24 \\
	\hline
	MUSYC1255:&  & &\\
	3.6 &  22.60 & 22.87 & 23.10  \\
	4.5 &  22.44 & 22.75 & 23.00   \\
	\hline
	\enddata
	\tablenotetext{a}{Using a 3$\arcsec$ circular aperture diameter.}
	\tablenotetext{b}{The 3$\sigma$ depths of this field  are taken from \citealt{ashby2018}.}
	\tablenotetext{c}{The  3$\sigma$ depths of this field are taken from \citealt{damen2011}. }

\end{deluxetable*}

\begin{figure*}[ht]
	\centering
	\includegraphics[scale=0.4]{./figures_journal/empty_ap_mag_cfhtlsd1.pdf}
	\caption{CFHTLSD1: Top row: 3$\sigma$ magnitude depths of $\mathrm{3.6\, \mu m}$ and $\mathrm{4.5\, \mu m}$ mosaics. Bottom row: Histograms of enclosed fluxes measured on $\sim$ 8000 empty apertures in $\mathrm{3.6\, \mu m}$ and $\mathrm{4.5\, \mu m}$ on the noise-normalized images. Best-fit Gaussians are also shown. }
	\label{f:empty_cfhtlsd1}
\end{figure*}

\begin{figure*}[ht]
	\centering
	\includegraphics[scale=0.4]{./figures_journal/empty_ap_mag_cfhtlsd3.pdf}
	\caption{ Same as in Figure~\ref{f:empty_cfhtlsd1}, but for the CFHTLS-D3 field.}
	\label{f:empty_cfhtlsd3}
\end{figure*}

\begin{figure*}[ht]
	\centering
	\includegraphics[scale=0.4]{./figures_journal/empty_ap_mag_cfhtlsd4.pdf}
	\caption{ Same as in Figure~\ref{f:empty_cfhtlsd1}, but for the CFHTLS-D4 field.}
	\label{f:empty_cfhtlsd4}
\end{figure*}

\begin{figure*}[ht]
	\centering
	\includegraphics[scale=0.4]{./figures_journal/empty_ap_mag_musyc1030.pdf}
	\caption{MUSYC1030: Top row: 3$\sigma$ magnitude depths of $\mathrm{3.6\, \mu m}$ and $\mathrm{4.5\, \mu m}$ mosaics. Bottom row: Histograms of enclosed fluxes measured on $\sim$ 4000 empty apertures in $\mathrm{3.6\, \mu m}$ and $\mathrm{4.5\, \mu m}$ on the noise-normalized images. Best-fit Gaussians are also shown. }
	\label{f:empty_musyc1030}
\end{figure*}

\begin{figure*}[ht]
	\centering
	\includegraphics[scale=0.4]{./figures_journal/empty_ap_mag_musyc1255.pdf}
	\caption{Same as in Figure~\ref{f:empty_musyc1030}, but for the MUSYC1255 field. }
	\label{f:empty_musyc1255}
\end{figure*}

\section{Public Data release}\label{s:release}
The data release consists of reduced images of all IRAC observations in three CFHTLSD fields, namely CFHTLS-D1, CFHTLS-D3, CFHTLS-D4,  and two MUSYC fields, MUSYC1030, and MUSYC1255. The images are available on the website: \url{http://cosmos.phy.tufts.edu/~danilo/NMBS2irac}.
The data release contains the following:
\begin{itemize}
\item Mosaic science images in both 3.6 $\mu$m and 4.5 $\mu$m. The images are W-registered to the reference image. This reference image is either the z-band stack image of the 25\% best seeing images from the CFHTLSDeep for CFHTLS-D1, CFHTLS-D3, CFHTLS-D4, or the K-band images obtained at the Blanco 4-m telescope at CTIO for MUSYC1030 and MUSYC1255.
\item Coverage maps containing exposure times (seconds) in both 3.6 $\mu$m and 4.5 $\mu$m.
\item Maps with 3$\sigma$ magnitude depth (as derived in Section~\ref{s:prop}) at the same WCS and pixel scales of the scientific images.
\end{itemize}

\section*{Acknowledgment}
The reduction and public release of the Spitzer-IRAC data in the NEWFIRM Medium-Band Survey II was supported by the National Aeronautics and Space Administration (NASA) under award number NNX16AN49G issued through the NNH15ZDA001N Astrophysics Data Analysis Program (ADAP). The Cosmic Dawn Center is funded by the Danish National Research Foundation. We thank the anonymous referee for his/her useful comments.

\bibliographystyle{aa}