\section{Introduction}
Recent successes of convolutional neural networks (CNNs) have led to the use of deep learning in real-world applications.
In order to manipulate these deep learning models, deep CNNs are trained using multi-class datasets to find manifolds separating different classes well.
To meet this need, deep and parameter-rich networks have emerged that have the power to find manifolds for {a large number} of classes.
However, these deep CNNs suffer from the problem of overfitting due to their great depth and complexity, which results in a drop of performance at the test time.
In fact, even a small ResNet applied for a dataset such as CIFAR-100 \cite{cifar100} will easily overfit
with the converged train losses, whereas the test accuracy is significantly {lower}.
These phenomena have led to the need for learning DNN models with appropriate regularization to allow them to generalize better.
Regularizing a model to achieve high performance for new inputs is a technique that has been used since the era of early machine learning.

\begin{figure}[t]

    \includegraphics[width=0.46\textwidth]{problem.png}
    \caption{Our problem formulation. The figure at the top shows how the KD is trained using ensemble of networks as the teacher.
    The figure at the bottom shows that a problem arises when we want to distill ensemble information directly in feature level rather than using label information.
    \label{fig:introduction}}

\end{figure}

Model ensemble \cite{dietterich2000ensemble} is one of the popular regularization methods, which  has been used as a way of alleviating the problem of overfitting in a single model.
However, it has drawbacks in that it requires multiple models and inputs should be fed to each of them at the test time.
Many studies proposed ideas to transfer knowledge of {a} teacher to {a} compact student~\cite{ba2014deep,rusu2015policy,modelcompression,urban2016deep}.
For a solution to this problem, \cite{hinton2015distilling} proposed {\textit{Knowledge Distillation}} (KD) which trains a student network using soft labels {from} an ensemble of multiple models or a {teacher network with high capacity}.
{They obtained} meaningful results in the speech recognition problem, and KD has become one of the representative methods of {\textit{knowledge transfer}}.
{They {aim} performance improvements of a weak student network by giving various forms of knowledge of expert teacher networks.}
It is also categorized as one family of model compression since it helps the student network achieve higher accuracy with a fixed number of parameters.

The recent knowledge transfer {algorithms} can be approximately categorized in two {viewpoints}.
The first is whether to use an ensemble model or a single high-capacity model as a teacher.
{The second is whether the teacher delivers predicted labels or the information from feature maps.}
Whereas the methods that use teacher's prediction can use both types of teachers, to the best of our knowledge, methods using {feature-map-level} information can only use a single high-capacity model as a teacher.
For example, studies of {\textit{Factor Transfer} (FT) \cite{kim2018paraphrasing}, \textit{Attention Transfer}  (AT)  \cite{zagoruyko2016paying}, and \textit{Neuron Selectivity Transfer}}  (NST) \cite{huang2017like} set the student network as {a} {shallow} {one} with a small number of parameters, and set the teacher network as {a} deeper and more powerful {one} instead of {an} ensemble expert.
One of the drawbacks of methods with a high-capacity teacher model is that high-capacity model may be hard to obtain \cite{lan2018knowledge}.

On the contrary, the ones which use {an} ensemble of networks can transfer ensemble knowledge and also have advantages of the peer-teaching framework \cite{furlanello2018born,hinton2015distilling,lan2018knowledge}.
Also, \cite{zhang2017deep} showed that using {the} same type of network {for transferring output-level knowledge} can improve the performance of a network.
{However, the methods that deliver knowledge at the feature-map level also have advantages that they can give more specific information to the student compared to the methods that only rely on the output predictions of {the} teacher.}

{To make full advantages of both the ensemble teacher method and the feature map transferring} method, we propose a new framework that {delivers} knowledge of multiple {networks} at {the} feature-map-level (See Fig. \ref{fig:introduction}).
In this paper, we train a new student network using multiple teachers that share the same architecture with the student network, named parallel FEED.
We also explore a {variant}, named sequential FEED, which recursively trains a new teacher.
Comparisons with other algorithms are provided with analysis that provides lessons.
Our main {contributions} are threefold:
{
\begin{itemize}
\item In various settings, we empirically show that for high-capacity networks, feature-map-based methods that give more specific knowledge are stronger than label-based methods which give more abstract knowledge.
\item We propose parallel FEED, a method that allows multiple teacher networks to be used for knowledge distillation at the feature-map level with non-linear transformations.
\item For qualitative analysis, {by utilizing the autoencoder reconstruction losses, we provide} hints for the {performance} enhancements {of our FEED.}
\end{itemize}
}

{The paper is organized as follows.}
First, we briefly explain the related works in the area of knowledge distillation.
Then our main proposed method is described
and a variant version implemented to compare with other methods is proposed.
Next, we verify our proposed methods with experiments.
{Experimental} results from our proposed training methods are compared with AT, KD, FT, BAN on CIFAR-100.
The ImageNet dataset is also used to check the feasibility of our method on a large-scale dataset.
Qualitative analysis is also provided which is followed by discussion and conclusion.

\section{Related Works}
Many researchers studied the ways to train models other than {using a} purely supervised loss.
In the early times of these studies, \textit{Model Compression} \cite{modelcompression} studied the ways to compress information from ensemble models in one network.
Ba \cite{ba2014deep} showed that shallow feed-forward nets can learn the complex functions previously learned by deep neural nets, by {minimizing} ${L}_{2}$ loss between logits.

More recently, Hinton \cite{hinton2015distilling} proposed \textit{Knowledge Distillation} {(KD)}, which uses softened softmax labels from teacher networks in training the student network {by minimizing the following loss}:
\begin{equation}
\small
\mathcal{L}_{KD}=(1-\alpha)\mathcal{L}_{CE}(y, \sigma(\textbf{s})) + \alpha T^{2}\mathcal{L}_{KL}\bigg(\sigma\big(\frac{\textbf{s}}{T}\big), \sigma\big(\frac{\textbf{t}}{T}\big)\bigg),
\label{eq:L_KD-likeFT}
\end{equation}
where the $\sigma(\cdot)$ is the softmax function, $T$ is a temperature value that controls the softened logit. $\alpha$ is hyper-parameter that controls the weight between two terms.
The vectors $\textbf{s}$ and $\textbf{t}$ are predicted output logits of the student network and the teacher network, respectively, and $y$ is the ground-truth label.
$\mathcal{L}_{CE}$ is the cross-entropy loss that is commonly used in classification problems, and $\mathcal{L}_{KL}$ is the Kullback-Leibler divergence loss.

This motivated researchers to develop many variants of it to various domains, and many researchers studied ways to better teach the student \cite{frosst2017distilling,pham2018efficient,radosavovic2018data,tan2018learning,tarvainen2017mean}.
We introduce some recent research flows that are potentially related to our proposed method.

\noindent \textbf{Peer teaching framework:}
\sw{Rather recently, many papers adapt peer-teaching framework that use the same kind of network for \nj{both the} teacher and \nj{the} student.} Geras \cite{geras2015blending} try to transfer knowledge between two networks that have almost the same number of parameters.
BAN \cite{furlanello2018born} shows that using the exact same architecture for the teacher and the student boosts the performance of the student network even without softening the labels.
{BAN uses a simpler loss term without softening both logits which KD does, and does not even assign weights to the two terms as follows:}
\begin{equation}
\mathcal{L}_{BAN}=
\mathcal{L}_{CE}\big(y, \sigma(\textbf{s})\big) +
\mathcal{L}_{KL}\big(\sigma(\textbf{t}), \sigma({\textbf{s}})\big).
\label{eq:BAN}
\end{equation}
They train the student recursively to enhance the performance further. The $n^{th}$ student network becomes the $(n+1)^{th}$ teacher network to train the next student network. The better teacher network will {teach the student network better}.

Also, studies such as DML \cite{zhang2017deep} and ONE \cite{lan2018knowledge} use the same kind of network {for on-line training of peer networks} with mutual KL-divergence losses.

\noindent \textbf{Feature-map-based methods for knowledge transfer:}
Contrary to the methods that try to use labels from the teacher network, there exist studies that distill useful information directly from feature maps in various forms.

AT \cite{zagoruyko2016paying} tried to transfer the attention map of the teacher network to {the} student network, and got meaningful results in knowledge transfer and transfer learning tasks. Their loss term is:
\begin{equation}
\small
\mathcal{L}_{AT}=\mathcal{L}_{CE}\big(y, \sigma(\textbf{s})\big) + \beta\sum^{L}_{l=1}\|\frac{f({A}^{t}_{l})}{\|f({A}^{t}_{l})\|_2}-\frac{f({A}^{s}_{l})}{\|f({A}^{s}_{l})\|_2} \|_2 ,
\label{eq:L_KD-likeFT}
\end{equation}
where $\beta$ is a hyperparameter that depends on the number of elements, and $l$ denotes the $l^{th}$ group \cite{yosinski2014transferable} in the network. $A^{t}_{l}$, $A^{s}_{l}$ are attention maps obtained from the teacher network and the student network, and $f(A) = (1/N)\sum^{N}_{n=1}a^{2}_{n}$, where $N$ is the number of channels, and $a_n$ is the spatial map from the $n^{th}$ channel.

Yim \cite{yim2017gift} {introduced another} knowledge transfer technique for faster optimization and applied it also to transfer learning. Shen\cite{shen2019amalgamating} tried to combine knowledge trained from different domains, and You \cite{You:2017:LMT:3097983.3098135} utilized ensemble of orderings of samples to teach the student network, which is very novel.
Huang \cite{huang2017like} tried to match the {features of the student and teacher networks by}
devising a loss term MMD (Maximum Mean Discrepancy).

FT \cite{kim2018paraphrasing} uses additional paraphraser and translator networks which help training the student network and got meaningful results.
Their loss terms are:
\begin{equation}
	\mathcal{L}_{rec} = \| x-P(x) \|^2_2, \quad {\text{(for the paraphraser)}}
    \label{eq:L_AE}
\end{equation}
\begin{equation}
\mathcal{L}_{student}=\mathcal{L}_{CE}(y, \sigma(\textbf{s})) + \beta\|\frac{F_T}{\|F_T\|_2}-\frac{F_S}{\|F_S\|_2} \|_1,
\label{eq:L_student}
\end{equation}
where $P(\cdot)$ is the autoencoder-based paraphraser network and $x$ is the input feature map for the paraphraser.
$F_T$ and $F_S$ are the output of the paraphraser and the translator, respectively.

\section{Proposed Training {Algorithms}}
In deep CNNs, due mainly to the curse of dimensionality, the data points that lie on the data space are very sparse. {For example, CIFAR datasets that contains the number of 50,000 training images and has 3,072 dimensions, so distances between each samples are very far.}
{Necessarily, decision boundaries that determine the borders dividing {classes} are multitudinous, because finding boundaries that {fit} well to {a training} dataset is {relatively} an easy task.
{Even if the {networks with the same architecture} are trained, the learned decision boundaries cannot be the same.}}
This is why ensemble methods usually {perform} better than a single model despite their structural equality.
{Goodfellow} \cite{Goodfellow-et-al-2016} also state that different models will not make all the same errors on the test dataset.

Consider the conditions that {determine} the training procedure of {CNNs}.
{They include} the structure of the CNN and the choice {of an} optimizer, the seed of random initialization, the sequence of mini-batches, and the {types of data} augmentations.
If one makes the same conditions for two different CNNs, their training procedure will be identical.
However, we usually determine only the structure of the CNN, usually keeping others to be random.
Consequently, two networks with {the} same structure will definitely not learn the same decision boundaries.

{Additionally,} {Kim \cite{kim2018paraphrasing}} {state} that they resolve the `inherent difference' between two networks.
Among the inherent differences, minimizing the differences in the structure of CNNs can help better learn the knowledge of the teacher network.
This has a thread of connection with that of {BAN \cite{furlanello2018born}}, which {shows} that using the same architecture for both {the} student and {the} teacher is actually beneficial.
This motivation provides us chances to produce several modified versions of existing methods.
In this section, we explain the {feature-level ensemble} training algorithms that are used for boosting the performance of {a} student network without introducing any {additional} calculations at {the} test time.
The proposed method is named as FEED which is an abbreviation for the \textit{FEature-level Ensemble for knowledge Distillation}.
We propose pFEED (\textit{parallel FEED}), which we handle as our main method, that use non-linear transformation layers to distill ensemble knowledge into the student network.
We also {introduce} sFEED (\textit{sequential FEED}), motivated by BAN \cite{furlanello2018born}, which adapts the sequential training with the use of nonlinear transformations.
The use of the nonlinear layers, rather than {using} a simple distance metric, had been explored previously in FitNet \cite{romero2014fitnets} and FT \cite{kim2018paraphrasing}.

\begin{figure}[t]
\centering
\includegraphics[width=1\linewidth]{pFEED2.png}\\
\caption{Illustration of our proposed {method}, parallel FEED. {NTL is an} abbreviation of Nonlinear Transformation Layer, and {one NTL is} allocated to each teacher {network}. All teacher networks are fixed during the training, and the student network and NTL networks are {trained simultaneously}.}
\label{pfeed}
\end{figure}

\subsection{{Parallel FEED}}\label{sec:pfeed}

Assuming that {distilling} knowledge {at the} feature map level and {distilling} knowledge {from an ensemble of multiple networks} both have their advantages, we wanted to make cooperation of these two kinds of methods.
Tackling this problem, we propose a training structure named \textit{pFEED} to transfer the ensemble knowledge at the feature map level.
Let us denote the number of the teacher network as ${N}$.
We use {different} non-linear transformation layers {for} each {teacher}, which means if there are ${N}$ teachers,
there are ${N}$ different nonlinear layers.
The final feature map of the student network is fed into the non-linear layers and its output is trained to mimic the final feature map of the teacher network.
In this way, we take {advantages} of both the ensemble model and the feature-based method.
Our {method} is illustrated in Figure \ref{pfeed}.
If we use ${N}$ different teachers, the loss term {is as follows:}
\begin{equation}
\mathcal{L}_{student}=\mathcal{L}_{CE}(y, \sigma(\textbf{s})) + \beta\sum^{N}_{n=1}{\mathcal{L}_{FEED_{n}}},
\label{eq:L_KD-likeFT}
\end{equation}
\begin{equation}
\mathcal{L}_{FEED_{n}} = \|{\frac{x^T_{n}}{\|x^T_{n}\|_2}-\frac{{NTL_n}(x^S)}{\|{NTL_n}(x^S)\|_2}}\|_1.
\label{eq:L_FT}
\end{equation}
{Here, $\mathcal{L}_{CE}$ is the} cross-entropy loss, $y$ {is the} ground-truth label, $\textbf{s}$ {is} the predicted logits, {and $\sigma(\cdot)$ denotes the softmax function.} $\mathcal{L}_{FEED_{n}}$ is the FEED loss from {the} $n^{th}$ teacher network, $x^T_n$ is the output feature map obtained from {the} $n^{th}$ teacher network, and ${x^S}$ is the output feature map obtained from the student network. ${NTL_n}(\cdot)$ is {the} $n^{th}$ nonlinear transformation layer used for {adapting the student with the} $n^{th}$ teacher network. \sw{Each ${NTL_n}(\cdot)$ is composed of three convolution layers with \nj{the} kernel size of 3 to expand the size of receptive field, so that the student can flexibly merge the knowledge attained from different teachers.}
The feature maps are normalized by its own size as in (\ref{eq:L_FT}). This normalizing term was previously used in AT. {$\beta$} is used to scale the $\mathcal{L}_{1}$ distance loss to match the scale of $\mathcal{L}_{CE}$, {as also described in AT}

\begin{figure}[t]
\centering
\includegraphics[width=1\linewidth]{sFEED3.png}\\
\caption{The sequential FEED. We adapted the setting of BANs that the teacher network and the student network have the same architecture
 The trained student network is used as a teacher network for the next {stage}.
{Dots mean} that we repeat this procedure.
}
\label{fig:sfeed}
\end{figure}

\subsection{{Sequential FEED}}\label{sec:sfeed}
BAN \cite{furlanello2018born} used cross-entropy loss combined with KD loss without softening the softmax logits.
They use a trained student network as a new teacher {which is used to train a new student} and do this recursively.
{We take this architectural advantage of using the same type of network recursively}
because it {is a suitable model for accumulating and assembling knowledge.}
{By} performing knowledge transfer several times recursively, it may also ensemble knowledge of many training sequences.
We applied the FEED recursively and named this framework as \textit{sFEED (sequential FEED)}.
{The training procedure of sFEED} is illustrated in Figure \ref{fig:sfeed}.

If the student network is trained standalone, it would perform similarly to the teacher network.
However, from the view of knowledge ensemble, since {the} teacher network delivers {feature-level} knowledge different from {that of the} student network, the student network will benefit from it.

\section{Experiments}
\label{sec:exp}
Firstly, we want to briefly show with experiments that using non-linear transformations for the output feature map with a distance metric is helpful for teaching the student network, rather than just {using} a simple distance metric {without any adaptation layer}.
Next, we will report the score of pFEED, which use multiple pairs of nonlinear transformation layers with each teacher networks, and will compare ours with KD and AT in {a} similar setting.
Finally, sFEED will be compared with FT, KD, and BAN.
The algorithms will be experimented following the BAN's sequential training schedules.

We show the classification results on CIFAR-100 \cite{cifar100} {on} which many networks show lower test accuracy than train accuracy that many studies do experiments on it to show their regularization power.
On this dataset, our results are compared with feature map based methods and label based methods with {corresponding} settings.
Second, we explore the feasibility of our algorithm on Imagenet \cite{ILSVRC15}, a commonly used large dataset and analyze the results quantitatively. In the remaining section, we show some analysis of our algorithm. The implementation details are on the supplementary material.

We chose three types of CNNs to check the applicability of our algorithms on CIFAR-100: ResNet \cite{he2016deep}, Wide ResNet \cite{zagoruyko2016wide}, and ResNext \cite{xie2017aggregated}.
For ResNets, we chose ResNet-56 and ResNet-110 which have fewer number of parameters compared to recent CNNs. WRN28-10 is a model that controls the widen factor, with much more number of parameters.
WRN28-10 achieves the best classification accuracy on CIFAR-100 among the WRNs reported {in \cite{zagoruyko2016wide}}.
The ResNext29-16x64d also achieves the best classification accuracy on CIFAR-100 in {\cite{xie2017aggregated}.
This} type of CNNs controls the cardinality of CNNs, and it has much more parameters compared to other models.
For ImageNet, we used ResNet-34 to confirm the feasibility {of our method} on large scale {datasets}.

\subsection{{Effectiveness of FEED loss}}
We compare the classification results using our FEED loss and two other {kinds} of loss terms, AT and simple $\mathcal{L}_{1}$ loss, {in training} the student at {the feature map level}.
The results are shown {in} Table \ref{TOY}.
The $\mathcal{L}_{1}$ means we simply use $\mathcal{L}_{1}$ loss at the final feature maps.
AT use attention maps attained from feature maps to give information.
From these {results}, we can see that using nonlinear transformation layers are helpful {in} delivering information at the feature-map level.
Interestingly, AT and $\mathcal{L}_{1}$ {beats} {a} single FEED model in ResNext29-16x64d model.

\begin{table}[t]
\caption{
Test classification error on {Cifar-100} dataset.
The numbers on the \textbf{scratch} column {are the} baseline errors of each network, trained by the pure cross-entropy loss, which {are the scores} of the teachers. The numbers on \textbf{scratch*} columns are the reported errors on their original papers.
}
\resizebox{1.\linewidth}{!}{
{\renewcommand{\arraystretch}{1.2}
\begin{tabular}{l|cc|ccc}
\Xhline{3\arrayrulewidth}
\multicolumn{1}{c|}{Model Type} & scratch*& scratch & $\mathcal{L}_{1}$  & AT    & FEED  \\%& pAT   & pFEED \\
\hline
ResNet-56                       & -& 28.18   & 27.16              & 26.60 & 26.02 \\%& 27.01 & \textbf{24.64}  \\
ResNet-110                      & -& 26.97   & 25.42              & 25.70 & 25.25 \\%& 25.23 & \textbf{22.84}  \\
WRN28-10                        & 19.25& 19.09   & 17.94              & 17.86 & 17.68 \\%& 17.73 & \textbf{17.10}  \\
ResNext29-16x64d                & 17.31& 17.32   & 16.46                  & 16.51 & 16.80 \\%& -     & \textbf{15.62}  \\
\Xhline{3\arrayrulewidth}
\end{tabular}
}
}
\label{TOY}
\end{table}

\subsection{Parallel FEED}\label{ex:parallel}
In {this experiment} on {pFEED}, our main experiment, we used the same type of networks that were used on previous experiments.
For all four types of CNNs, we compared the classification results with the result of KD because we designed our training algorithm with the intention of {distilling} more ensemble-like knowledge from multiple teachers.
{We also experimented AT to empirically show that it would not be easy for a feature-map-based knowledge transfer methods to fully utilize multiple teacher networks}.
We modulated AT to use multiple teachers, by simply incrementally adding same loss terms for each teacher network. We did not change the $\beta$ values for weights of each loss terms, {similar} to pFEED.

The results are in Table \ref{pfeed}.
The `Scratch' column shows the performance of the base networks, {which are} used in KD for model ensemble, and also used as teachers in pFEED and AT.
For all experiments, pFEED consistently got higher accuracy compared with both KD and AT and produces the closest {results} to {those of the} network ensemble.

\begin{table*}[]
\caption{{Test classification error (\%) on CIFAR-100 dataset. All scores of other methods are our reproduction. In the 5x$\mathcal{L}_{1}$, 5xAT and 5xFT columns, we used five $\mathcal{L}_{1}$, AT and FT losses each for training one student.}
The {scores of the Ens} column {are the} performance of label ensemble of 5 scratch models. We trained pFEED 5 times and averaged the results.}
\begin{center}
\begin{tabular}{l|c|ccccccc}
\Xhline{3\arrayrulewidth}
\multicolumn{1}{c|}{Model Type}           &\multicolumn{1}{c|}{Scratch (5 mean)}   & KD & 5x$\mathcal{L}_{1}$    &  5xAT   & 5xFT & pFEED  &  Ens & Parameters \\
\hline
ResNet-56            & 28.18 ($\pm$ 0.17) & \textbf{24.69} & 27.29 & 27.02 & 25.41 & 24.74 ($\pm$ 0.12) & 22.45 & 0.85M\\
ResNet-110           & 26.97 ($\pm$ 0.16) & 23.50 & 28.94 & 25.23 & 23.81 & \textbf{22.98} ($\pm$ 0.20) & 21.20 & 1.73M\\
WRN28-10             & 19.09 ($\pm$ 0.13) & 18.30 & 23.75 & 17.73 & 17.30 & \textbf{16.86} ($\pm$ 0.14) & 16.59 & 36.5M\\
ResNext29-16x64d     & 17.32 ($\pm$ 0.08) & 16.64 & 16.31 & 17.42 & 16.30 & \textbf{15.70} ($\pm$ 0.08) & 15.66 & 68.1M\\
\Xhline{3\arrayrulewidth}
\end{tabular}
\end{center}
\label{pfeed}
\end{table*}

\textbf{Comparison with KD:} It is worth noting that the performance of KD is almost equivalent to pFEED for small networks, but as it comes to the networks with {a} larger number of parameters, pFEED shows better accuracy compared to KD.
This result matches the hypothesis that {distilling in the feature-map level} will {provide} more \textbf{detailed} information {to the student}.
{In contrast, {KD} works quite well for small networks, because {it gives} ensemble labels, which {are} rather \textbf{abstract}.
These labels are useful for small networks that should focus on key information {for accuracy improvement}.}

\sw{\textbf{Comparison with feature-map-based methods:}  Compare the results on Table \ref{pfeed} with Table \ref{TOY}. As shown by the scores of $\mathcal{L}_{1}$ and AT, using multiple teachers \nj{with these methods} did not make a meaningful difference compared to the case when a single teacher network was used.}
{This is} {an example} of our statement, that it is not easy for existing {feature-map-based} knowledge transfer methods to {fully} utilize multiple teachers to boost the performance of a single student network.

The results of pFEED for ImageNet is on Table \ref{table:pfeed-imagenet}.
We also could find some accuracy improvements on ImageNet dataset, but did not have enough resources to train models with larger parameters, nor did we could experiment on larger models or other methods.
Though we reported scores of five scratch models, only {the} first three teacher networks had been used for training {of the student}.
We could get decent results, but interestingly, improvements are not as strong as those of sFEED on ImageNet that will be shown later.

\begin{table}[t]
\caption{Validation classification error (\%) of pFEED on Imagenet dataset.}
\resizebox{1.\linewidth}{!}{
{\renewcommand{\arraystretch}{1.2}
\begin{tabular}{l|ccccc|c}
\Xhline{3\arrayrulewidth}
\multicolumn{1}{c|}{Model Type}  &\multicolumn{5}{c|}{Scratch (5 runs)} & pFEED \\
\hline
ResNet-34(Top-1)   & 26.45    & 26.59 & 26.40 & 26.77 & 26.64 & \textbf{25.27}  \\
ResNet-34(Top-5)   & 8.54     &  8.72 &  8.63 &  8.68 &  8.61 & \textbf{7.79}  \\
\Xhline{3\arrayrulewidth}
\end{tabular}
}}
\label{table:pfeed-imagenet}
\end{table}

\subsection{Sequential FEED}\label{sec:compare}
This section contains results of 4 types of {methods (FT, KD, BAN, and sFEED)} experimented in small {and large} capacity networks.
Note that the pFEED is our main proposed method for good performance, and {the} purpose of {this subsection} is to further delve into each method to speculate the characteristics and differences of each {method}.
The classification results of {different methods} on CIFAR-100 can be found in Table \ref{table:sfeed_comp}.
The word `Stack' in {Table} \ref{table:sfeed_comp} is the number of recursions that the student model is trained.
We only experimented up to 5 times since all of them {achieve} fairly good enough accuracy compared to baseline models.
We reported the {performance of} sFEED because it has a more similar framework to BANs, though the pFEED, our main proposed method, outperforms {sFEED}.
The BAN-N in {\cite{furlanello2018born}} would be {the} identical setting to the Stack-(N-1) in Table \ref{table:sfeed_comp}.
We additionally experimented FT with {the identically structured teacher and student networks}.
FT basically uses {a} large teacher network {with} a paraphraser, which is trained as an autoencoder.
It is trained to extract key information called `factor' from the teacher network in {an} unsupervised {manner}, so it gives more abstract information.

\begin{table*}[]
\caption{Test classification error (\%) of sFEED on CIFAR-100 dataset. The model's scores on the \textbf{Scratch*} column are the same as the scores reported on their original papers, and {those} on the \textbf{Scratch} column are from our implementation. The parameters are counted in Millions.}
\begin{center}
\resizebox{0.7\linewidth}{!}{
\begin{tabular}{c|l|c|c|c|ccccc}
\Xhline{3\arrayrulewidth}
\multirow{16}*{CIFAR-100}& \multicolumn{1}{c|}{Model Type}  & Scratch             & Distillation Type & Algorithm & Stack-2& Stack-3& Stack-4& Stack-5&  \\
\cline{2-10}
                        & \multirow{4}*{ResNet-56}          & \multirow{4}*{28.03}& \multirow{2}*{Label}   & BAN        &25.85 & 25.52 & \textbf{25.30} & 25.45 & \\
                        &                                   &                     &                        & KD         &25.55 & 25.14 & \textbf{24.98} & 24.98 & \\
\cline{4-5}
                        &                                   &                     & \multirow{2}*{Feature} & FT         &26.33 & 25.70 & 25.91 & \textbf{25.18} & \\
                        &                                   &                     &                        & sFEED      &26.02 & 26.00 & 25.59 & \textbf{25.33} & \\
\cline{2-10}
                        & \multirow{4}*{ResNet-110}         & \multirow{4}*{27.14}& \multirow{2}*{Label}   & BAN        &24.67 & \textbf{23.81} & 23.98 & 24.06 & \\
                        &                                   &                     &                        & KD         &24.35 & 24.01 & 23.73 & \textbf{23.66} & \\
\cline{4-5}
                        &                                   &                     & \multirow{2}*{Feature} & FT         &24.90 & 24.50 & 24.34 & \textbf{24.13} & \\
                        &                                   &                     &                        & sFEED      &25.25 & \textbf{24.33} & 24.58 & 24.40 & \\
\cline{2-10}
                        & \multirow{3}*{WRN28-10}           & \multirow{3}*{19.00}& \multirow{1}*{Label}   & KD         &18.47 & \textbf{18.43} & 18.57 & 19.05 & \\
\cline{4-5}
                        &                                   &                     & \multirow{2}*{Feature} & FT         &18.23 & 18.05 & 18.14 & \textbf{17.84} & \\
                        &                                   &                     &                        & sFEED      &17.68 & 17.50 & 17.52 & \textbf{17.27} & \\
\cline{2-10}
                        & \multirow{4}*{ResNext29-16x64d}   & \multirow{4}*{17.31}& \multirow{2}*{Label}   & BAN        &16.59 & \textbf{16.42} & 16.43 & 16.48 & \\
                        &                                   &                     &                        & KD         &\textbf{16.87} & 16.90 & 16.88 & - & \\
\cline{4-5}
                        &                                   &                     & \multirow{2}*{Feature} & FT         &16.80 & 16.76 & \textbf{16.47} & 16.48 & \\
                        &                                   &                     &                        & sFEED      &16.80 & 16.47 & 16.22 & \textbf{15.94} & \\
\Xhline{3\arrayrulewidth}
\end{tabular}
}
\end{center}
\label{table:sfeed_comp}
\end{table*}

{\textbf{Comparison of Label-based-methods and Feature-map-based methods:}
For the smaller networks such as ResNet-56 and ResNet 110, methods using labels performed better than the {feature-map-based} methods, but for networks with larger sizes, {feature-map-based} methods showed higher accuracy.
FT {uses} more abstract knowledge (using paraphraser) compared to sFEED, so it did not perform well as sFEED for larger networks. {However,} it performed better for smaller-sized networks.
KD use{s} more abstract knowledge compared to BAN, because KD soften{s} the labels, and it achieve{s} higher accuracy for smaller networks, but BAN showed better accuracy for ResNext.}

The results of {sFEED} for ImageNet is on Table \ref{ex-imagenet}.
For the base model, we simply used the pre-trained model that Pytorch supplies, and could achieve {the desired} result that the performance of Top-1 and Top-5 accuracy improves at each {Stack}.
The sFEED with Stack-5 achieves better performance compared to pFEED. {The performance of} pFEED in Table \ref{table:pfeed-imagenet} trained with only three teachers {is} close to {that of} sFEED with Stack-3, which is a reasonable comparison.

\begin{table}[]
\caption{Validation classification error (\%) of sFEED on Imagenet dataset. The model's scores on the \textbf{Scratch*} column are the same as the scores reported on the Pytorch implementation.}
\resizebox{1.\linewidth}{!}{
{\renewcommand{\arraystretch}{1.2}
\begin{tabular}{l|cccccc}
\Xhline{3\arrayrulewidth}
\multicolumn{1}{c|}{Model Type}  &\multicolumn{1}{c}{Scratch*} & Stack-2 & Stack-3 & Stack-4 & Stack-5
\\ \hline
ResNet-34(Top-1)   & 26.45    & 25.60 & 25.30 & 25.18 & \textbf{25.00}  \\
ResNet-34(Top-5)   & 8.54     & 8.08  &  7.86 &  \textbf{7.73} &  7.83  \\
\Xhline{3\arrayrulewidth}
\end{tabular}
}}
\label{ex-imagenet}
\end{table}

\subsection{Qualitative Analysis}\label{recon}
\textbf{Reconstruction Loss:}
{Suppose} that the reason for the accuracy gains shown in the previous tables is that the student learns the ensemble knowledge
{that contains information with high complexity.}
But how can one actually distinguish whether the networks learn complex information or not?
Here, we adopt the paraphraser from FT.
A paraphraser in the FT is a convolutional autoencoder in that it uses convolution {and} transposed convolution layers {with a} reconstruction loss. Then the factor can be interpreted as a latent vector ${z}$.

Let us denote the input of paraphraser as {$x$}. The increase in the complexity of feature representation is equivalent to the increase in the complexity of $x$.
Since the number of parameters in {the} paraphraser is fixed, the size of $z$ also should be fixed.
Consequently, {as the complexity of $x$ increases}, $p(x|z)$ decreases, resulting in an increase of {the} reconstruction loss.
Reconstruction errors were used as a criterion for feature selection or PCAs in  \cite{boutsidis2008unsupervised,farahat2011efficient,li2017reconstruction,masaeli2010convex}, where they use linear models.
In our experiment, we use an arbitrary autoencoder composed of convolution layers with nonlinear activation.
\begin{figure}[t]
\centering

    {\includegraphics[width=0.3\textwidth]{20190813_sfeed.png}
    \caption{Paraphraser reconstruction loss $\mathcal{L}_{rec}$(training) for ResNet-56 with sFEED.}\label{rec_loss_step}}
    {\includegraphics[width=0.3\textwidth]{20190813_pfeed.png}
    \caption{Paraphraser reconstruction loss $\mathcal{L}_{rec}$(training) for ResNet-56 with pFEED.}\label{rec_loss_kd}}

\end{figure}

For sFEED and pFEED, we recorded the {average training} reconstruction losses of the paraphrasers {for ResNet-56} {normalized by the size of the paraphraser} and plotted the curve on Figure \ref{rec_loss_step} and Figure \ref{rec_loss_kd} (though we do not actually use paraphrasers for FEED training).
{In Fig. \ref{rec_loss_step}, St1} through St4 are the paraphrasers trained based on the student networks {of} Scratch {(St1)} through Stack-4 {(St4)} in Table \ref{table:sfeed_comp}.
The paraphrasers in Figure \ref{rec_loss_kd} is trained based on one of the scratch teacher networks and following student network {of pFEED {and KD}} in Table \ref{pfeed}.
{As expected, as the knowledge is transferred, the reconstruction loss becomes larger which indicates that the student network learns more difficult knowledge and thus the classifier accuracy increases.}
{This} trend matches the results on the tables.
In Figure \ref{rec_loss_step}, the {the trend of reconstruction losses reciprocally matches the accuracy trend} of ResNet-56 row in Table \ref{table:sfeed_comp} of sFEED.
{(Especially, Stack 2 and 3 have similar errors, and likewise, St2 and St3 are similar)}.
The big {difference in reconstruction loss between `Scratch' and `pFEED'} in Figure \ref{rec_loss_kd} also {corresponds to} the high performance increase in ResNet-56 row of Table \ref{pfeed}.
A better teacher network would learn more complex and detailed features because it {has} to contain powers to distinguish important {but} different details from each image to form better decision boundaries.
{It is worth noting that the curves for KD and pFEED shows opposite aspect, even though they both succeed in enhancing their performance.
Here, our assumptions also holds:
First, training the student with multiple teacher's feature map will help the student learn detailed features.
Second, teacher's labels are abstract information, but will help the student learn key information.}
Thus, the tendency of KD in Figure \ref{rec_loss_kd} is opposite to pFEED.
{The curves for other type of networks also show consistent aspect, and more examples are handled in the {supplementary} materials.}

\section{Discussion}
Experiment in Section \ref{ex:parallel} with pFEED shows that allocating nonlinear transformations for each of teacher networks can extract ensemble knowledge from multiple teachers.
{On the other hand,} AT, which directly mimics the attention map, fails to learn from multiple teacher networks.
The error of pFEED comes closer to {the} actual model ensemble compared to KD, especially for the models with high capacity.
Next, experiments in Section \ref{sec:compare} compares sFEED with FT, KD, and BAN. Results give lessons to the choice of the algorithm that would be useful depending on the type of networks.
The FT which extract key information from the teacher network performs better than sequential FEED for smaller networks and worse for larger networks.
The KD and BANs, {using} labels which is even more abstract, {perform} better than sequential FEED for smaller networks.
However, the result shows that sFEED with nonlinear transformation layers are more useful for networks {with a} larger capacity.
Though not absolute, if one wants to use distillation for model compression with smaller networks, it would be beneficial to seek ones that use abstract information like labels.
If one wants to use distillation for high performance where higher performance is needed, distillation methods that can give more detailed information can be useful.

The analysis of the reconstruction error, which we utilized the paraphraser of FT, would be helpful to judge whether the network compactly learned its features.
Finally, since we did not actively search better options for the nonlinear transformations, better choices for them might exist.

\section{Conclusion}
In this work, we proposed a couple of new network training algorithms referred to as \textit{FEature-level Ensemble {for} knowledge Distillation} (FEED).
With FEEDs, we can improve the performance of a network by trying to inject ensemble knowledge {in the feature-map level} to the student network.
The first one, parallel FEED trains the student network using multiple teachers simultaneously.
The second one, sequential FEED recursively trains the student network and incrementally improves performance.
The qualitative analysis with reconstruction loss gives hints about the cause of accuracy gains.
{The main drawback is the training times needed for multiple teachers which is an inherent characteristics of any ensemble methods,} and pFEED causes bottleneck by feeding inputs to multiple teachers simultaneously.
But, it does not affect the inference, which is beneficial without trade-offs in {the} test time.
{Devising a more train-efficient method} will be our future work, together with an application to other domains {other than classification tasks}.

\bibliographystyle{aaai}


\section{Introduction}

Intensity mapping of redshifted emission lines \citep{bharadwaj2001, loeb2008}
is a novel technique that has
the potential to perform precision cosmology by detecting the integrated
emission from sources across redshifts without resolving individual systems
\citep[e.g.,][]{kovetz2019}. Besides offering rich insights into the physics of
star formation history and the processes governing galaxy evolution
\citep[e.g.,][]{wyithe2008, wolz2016}, it has the
ability to improve vastly upon the current measurements of cosmological
parameters \citep[e.g.,][]{bull2014}, as well as place competitive constraints
on inflationary scenarios and physics beyond the standard model
\citep[e.g.,][]{hall2013, camera2013, pourtsidou2016, masui2010}. The most well-studied
example of line-intensity mapping involves that of the redshifted 21-cm
emission of neutral hydrogen (hereafter, HI) which arises primarily in
star-forming galaxies and the intergalactic medium at low to moderate
redshifts.

Using line-intensity mapping techniques in synergy with other, more traditional
and established tracers of large-scale structure is crucial to unlock the true
potential of these surveys. It is known that  \citep[e.g.,][]{seljak2009, fonseca2015, fonseca2017, camera2015}
cross-correlations of several individual tracers of the cosmological structure
often offer several significant advantages over individual surveys. { The
systematic survey-specific effects are mitigated to a large extent, the noise
in the surveys is reduced, and the foregrounds and contaminants of individual
surveys are,
in most cases, uncorrelated and hence do not bias the cross-correlation
measurement. \footnote{ However, we note that the presence of foregrounds may greatly increase the variance, which is an important effect, and hence efficient foreground cleaning or avoidance techniques are still necessary in order to isolate the true signal.}} The cosmic variance can be mitigated in the measurement of
some of the cosmological parameters \citep[e.g.,][]{mcdonald2009, abramo2013}.

To this end, cross-correlating 21 cm intensity mapping surveys with optical
galaxies offers rich possibilities into exploiting the complementarity of both
approaches. The first intensity mapping detection of the redshifted 21 cm
emission at $z \sim 0.53-1.12$ with the Green Bank Telescope (GBT) was made in
cross-correlation with the DEEP2 optical galaxy survey \citep{chang2013}, and
has been followed up since then resulting an updated cross-power spectrum using
the WiggleZ survey at $z \sim 0.8$ \citep{masui13} and an upper limit on the
auto-power spectrum \citep{switzer13}, which was used to place the first
intensity mapping constraints on the product of the neutral hydrogen density
and bias parameter. Similarly, the cross-power spectrum between 2dF galaxies in
the southern hemisphere and the Parkes HI intensity field at $z \sim 0$ has
also been measured recently \citep{anderson2018}, offering insights into the
clustering of low redshift HI systems.

The auto-correlation power spectrum signal of high-redshift 21 cm in emission has yet to
be observed, although there are several experiments planned or in the final
stages of commissioning to achieve this goal. These include (i) the Canadian
Hydrogen Intensity Mapping Experiment
(CHIME)\footnote{https://chime-experiment.ca/}, (ii) the Hydrogen Intensity and
Real-time Analysis eXperiment
(HIRAX)\footnote{https://www.acru.ukzn.ac.za/~hirax/}, (iii) the BAO In Neutral Gas
Observations \citep[BINGO;][]{battye2012}, (iv) the TianLai experiment
\citep{chen2012} and the Five hundred metre Aperture Spherical Telescope
\citep[FAST;][]{smoot2017},
(v) the Meer-Karoo Array Telescope \citep[MeerKAT;][]{jonas2009} and
(vi) the Square Kilometre Array (SKA) Phase I MID
\footnote{http://www.ska.ac.za/}.

The optical surveys of key interest for cross-correlations with the 21 cm
surveys planned above include those with (i) the completed Dark Energy Survey
(DES)\footnote{https://www.darkenergysurvey.org/}, a photometric
galaxy survey over $z\sim 0.5 -1.4$,  cataloging hundreds of millions of
galaxies in the southern hemisphere,  (ii) the forthcoming Dark Energy
Spectroscopic Instrument (DESI), a spectroscopic survey which will target a few
tens of millions of galaxies in the northern sky over the redshift range $z
\sim
0-3$, measuring cosmological parameters and the growth of structure through
redshift space distortions, and (iii) future galaxy surveys conducted with the
Large Synoptic Survey Telescope (LSST)\footnote{www.lsst.org}, and the
space-based Euclid\footnote{www.euclid-ec.org} spectroscopic survey.

In \citet{pourtsidou2016}, synergies between a MeerKAT HI intensity mapping
survey and photometric galaxies from the  Dark Energy Survey (DES) have been
explored. Recent studies \citep{mona2019, amadeus2019} have illustrated the
ability of HI intensity mapping (with SKA and HIRAX) in cross-correlation with
photometric galaxy surveys (such as DES and LSST) to measure the gravitational
lensing magnification. In \citet{carucci2017}, the synergies between 21 cm SKA
I MID  and the Baryon Oscillation Spectroscopic Survey (BOSS)-like Lyman-alpha
surveys have been presented, which can constrain the bias of astrophysical
systems. In \citet{cosmicvisions2018}, various prospects for cross-correlating
21 cm intensity mapping and optical surveys have been explored, including with
QSOs observed by the DESI survey.  It has been shown \citep{chen2018} that a
combining a CMB Stage 4-like survey with 21 cm intensity mapping observations
from a SKA I MID like survey cross-correlated with DESI quasars can enable
precise measurements of the growth factor, and test the predictions of general
relativity on the largest scales.
\citet{witzemann2018} explore how synergies between a SKA I MID like survey and
a photometric LSST survey can mitigate the effects of cosmic variance, enabling
measurements of the bias ratio at large scales up to $\ell \sim 3$, and
\citet{hall2017} illustrate how peculiar velocity effects can be constrained
using the dipole of the redshift space cross-correlation between 21 cm and
optical surveys conducted with various experiments. In \citet{ballardini2019a}, the constraints for local primordial non-Gaussianity
has been studied with a SKA-like intensity mapping survey in
cross-correlation with photometric galaxy surveys (Euclid and LSST)
and CMB lensing.

In this paper, we build upon our previous forecasting analyses in
\citet[][hereafter Paper I]{hparaa2018} which focussed on auto-correlation 21
cm power spectra and extend these to the case of measurement of the
astrophysical and cosmological parameters using a combination of 21 cm and
optical galaxy surveys. We use the uncertainties in the parameters coming from
the combination of current measurements to set realistic priors on the HI
astrophysics. We consider three sets of surveys in this work: (i) a CHIME-like
survey overlapping with DESI in the northern hemisphere, (ii) a MeerKAT-like
survey overlapping with DES in the southern hemisphere and (iii) a SKAI MID
survey overlapping with the Large Synoptic Survey Telescope (LSST), again in
the southern hemisphere.

The paper is organized as follows. For modelling the HI distribution and
density profile, we use the halo model framework introduced in \citet{hpar2017}
and expanded upon in \citet[][hereafter Paper II]{hparaa2016}, which describes
the best fitting HI-halo mass relation and profile  constrained by the
currently available data. The bias and redshift distribution of the optical
galaxies are modelled
following
the treatment for the particular survey under consideration. These frameworks
are
briefly described in  Sec. \ref{sec:framework}. Using the cross-correlation
power spectrum thus derived, we compute the relative errors on the
astrophysical and cosmological parameters under consideration using a Fisher
forecasting formalism for the three survey sets in Sec. \ref{sec:surveys}.  We
comment on the comparison of these predictions to those from the corresponding
21 cm auto-correlation
constraints, and summarize our conclusions in Sec. \ref{sec:conclusions}.

\section{Framework for Fisher forecasts}
\label{sec:framework}
We use the halo model for neutral hydrogen (see Paper II) and build upon our
existing forecasts in Paper I, which had focused on the HI auto-correlation
surveys alone. The halo model framework consists of a prescription assigning
average HI mass to halo mass $M$ at redshift $z$, given by:
\begin{eqnarray}
M_{\rm HI} (M,z) &=& \alpha f_{\rm H,c} M \left(\frac{M}{10^{11} h^{-1}
M_{\odot}}\right)^{\beta} \nonumber \\
&\times& \exp\left[-\left(\frac{v_{c0}}{v_c(M,z)}\right)^3\right]
\end{eqnarray}
In the above formula, the free parameters are given by: (i) $\alpha$, the
average HI fraction relative to cosmic $f_{\rm H,c}$, (ii)   $\beta$, the logarithmic slope which represents the deviation from linearity of the
HI-halo mass prescription, and (iii) $v_{c0}$, which denotes the minimum virial
velocity below which haloes preferentially do not host HI.

To model the smaller scales in the HI power spectrum, we also need a
prescription for the profile of the HI as a function of radius, halo mass and
redshift, which is found to be well modelled by an exponential function { (Paper II, see also the observational results from, e.g., \citet{bigiel2012})}:
\begin{equation}
\rho(r,M) = \rho_0 \exp(-r/r_s)
\label{rhodefexp}
\end{equation}
with the scale radius $r_s$ given by:
\begin{equation}
 r_s = R_v(M,z)/c_{\rm HI} (M,z)
 \end{equation}
with $R_v$ being the halo virial radius and $c_{\rm HI}$ being the
concentration parameter of the HI systems, which is analogous to the
corresponding expression for dark matter:
\begin{equation}
 c_{\rm HI}(M,z) =  c_{\rm HI, 0} \left(\frac{M}{10^{11} M_{\odot}}
\right)^{-0.109} \frac{4}{(1+z)^{\gamma}}.
\end{equation}
The HI profile thus introduces two more free parameters through the
concentration parameter and its evolution: (i) $c_{\rm HI, 0}$ representing the
overall normalization and $\gamma$ which encodes the evolution of the function
with redshift.
In order to compute the nonlinear HI power spectrum, we need the Fourier
transform of the profile function, given by:
\begin{equation}
u_{\rm HI}(k|M) = \frac{4 \pi}{M_{\rm HI} (M)} \int_0^{R_v} \rho_{\rm HI}(r)
\frac{\sin kr}{kr} r^2 \ dr
\end{equation}
which allows us to write the power spectrum for HI intensity fluctuations as
the sum of the 1- and 2-halo terms:
\begin{equation}
P_{\rm HI} (k,z) = P_{\rm 1h, HI} + P_{\rm 2h, HI}
\end{equation}
with
\begin{equation}
P_{\rm 1h, HI} (k,z) = \frac{1}{\bar{\rho}_{\rm HI}^2} \int dM \  n(M) \ M_{\rm
HI}^2 \ |u_{\rm HI} (k|M)|^2
\label{power1h}
\end{equation}
and

\begin{eqnarray}
&& P_{\rm 2h, HI}  (k,z) =  P_{\rm lin} (k) \nonumber \\
&&\left[\frac{1}{\bar{\rho}_{\rm HI}} \int dM \  n(M) \ M_{\rm HI} (M) \ b (M)
\ |u_{\rm HI} (k|M)| \right]^2
\label{power2h}
\end{eqnarray}
with $P_{\rm lin}(k)$ being the linear matter power spectrum.

For  computing the power spectrum of the optical galaxies in the survey, we use
the expression:

\begin{equation}
P_{\rm gal} (k,z) = P_{\rm dm} (k,z) \ b^2_{\rm gal} (k,z)
\end{equation}

where $P_{\rm dm}$ is the dark matter power and the $b^2_{\rm gal} (k,z)$ denotes
the galaxy-galaxy bias factor. This factor changes according to the survey and
the type of galaxies under consideration.
This (scale-dependent bias) is modelled following the parameters given by
\citet{amendola2015} which is based on the $Q$-formula of \citet{cole2005}:
\begin{equation}
b_{\rm gal} = b_{\rm gal, ls} b_0 \left(\frac{1 + Q k^2}{1 + Ak} \right)^{1/2}
\end{equation}
where the values $b_0 = 1.3, A = 1.7, Q = 4.6$ are assumed not to vary with
redshift. The $b_{\rm gal, ls}$ term depends on the survey under consideration.

The angular power spectrum on the sky for HI is computed by using the standard
result:
\begin{eqnarray}
			 &&C_{\ell,\rm{HI}}(z,z') = \frac {2} {\pi} \int
d\tilde z\, W_{\rm HI}(\tilde z) \int d\tilde z'\, W'_{\rm HI}(\tilde z')
\nonumber \\
			&& \hskip-0.2in \times  \int k^2 dk\, \langle \delta_{\rm HI}
(\textbf{k}, z) \delta_{\rm HI} (\textbf{k}', z') \rangle
j_\ell(kR(\tilde z)) j_\ell(kR(\tilde z')),
	\label{eq:lintheory}
\end{eqnarray}
In the above expression, $\langle \delta_{\rm HI} (\textbf{k}, z) \delta_{\rm
HI} (\textbf{k}', z') \rangle$ is the ensemble average of the HI density
fluctuations at $(\textbf{k},z)$ and $(\textbf{k}', z')$ respectively. This is,
in general, not expressible purely in terms of the power spectrum of HI as
defined above, $P_{\rm HI} (k,z)$ evaluated at either of $\{z,z'\}$ since the
density field evolves with $z$. However, in many cases,
one can approximate this as $P_{\rm HI} (k,z_{m})$ where $z_m$ is the mean
redshift of the given bin. In what follows, we use $z$ and $z_m$
interchangeably.

In the above expression, the $W_{\rm HI}, W'_{\rm HI}$ are the HI window
functions at the redshifts
$z$ and $z'$, taken
to be uniform across the redshift bin considered, and $R(z)$ is the co-moving
distance to redshift $z$. We use a
top hat window function $W_{\rm HI}(z)$ with a width of $\Delta z  = 0.5$.

For a generic galaxy survey, calculation of the angular power spectrum yields
the expression:
\begin{eqnarray}
			 &&C_{\ell,\rm{gal}}(z,z') = \frac {2} {\pi} \int
d\tilde z\, W_{\rm g}(\tilde z) \int d\tilde z'\, W'_{\rm g}(\tilde z')
\nonumber \\
			& &\times  \int k^2 dk\, P_{\rm gal}(k,z)
j_\ell(kR(\tilde z)) j_\ell(kR(\tilde z')),
	\label{eq:lintheory}
\end{eqnarray}
The dark matter power spectrum for linear scales can alternatively be written
as $P_{\rm dm} (k,z)
= P_{\rm dm} (k,0) D^2(z)$ where
$D(z)$ is the growth factor for the dark matter
perturbations whose power spectrum is normalized such that $D(0) = 1$.
The window function for the galaxy survey, $W_{\rm g}$, can be different
from that of the HI, and depends on the details of the selection function
(usually denoted by $\phi(z)$) of each survey.
 Parametrized forms for $\phi(z)$ are available for different galaxy surveys
and usually follow a standard functional form \citep{smail1995}:
\begin{equation}
 \phi(z) \propto z^{\alpha} \exp(- (z/z_0)^{\beta})
 \label{selfunc}
\end{equation}
where $\alpha$, $\beta$ and $z_0$ are fitted from the galaxy counts data in
different redshift bins.
Once $\phi(z)$ is known, we derive the window function for the survey as:
\begin{equation}
 W_{\rm g}(z) = \phi(z)/\int_{z_{\rm min}}^{z_{\rm max}} \phi(z) dz
\end{equation}
where $z_{\rm min}$ and $z_{\rm max}$ are the redshift edges of the survey.
This ensures that the window function is normalized, i.e.
\begin{equation}
 \int_{z_{\rm min}}^{z_{\rm max}} W_g(z) dz = 1
\end{equation}

 The calculation of the angular power spectra above, both for HI and for
galaxies, can be simplified on using the
Limber approximation \citep{limber1953} which is a good approximation in the
large $\ell$ ($\ell > 50$) limit. The expression can be shown to reduce to:
\begin{equation}
C_{\ell, {\rm HI/gal}} = \frac{1}{c} \int dz  \frac{{W_{\rm HI/gal}}(z)^2
H(z)}{R(z)^2}
P_{\rm HI,gal} [\ell/R(z), z]
\label{cllimber}
\end{equation}.

The cross-correlation signal is then calculated as:
\begin{equation}
C_{\ell, \times} = \frac{1}{c} \int dz \frac{{W_{\rm HI}(z) W_{\rm gal}(z)}
H(z)}{R(z)^2}
(P_{\rm HI} P_{\rm gal})^{1/2}
\label{clcross}
\end{equation}
where the arguments of both power spectra ($P_{\rm HI}$ and $P_{\rm gal}$) are
at $[\ell/R(z), z]$.
{ Noise in the HI intensity mapping survey, for the MeerKAT-like and SKA I MID-like configurations, is calculated using the standard
expression assuming the interferometer array to operate in the single-dish
autocorrelation mode \citep[e.g.,][]{Knox1995,ballardini2019}:
\begin{equation}
N_{\ell, \rm HI} = \left(\frac{T_{\rm sys}}{\bar{T}(z)}\right)^2
\left(\frac{\lambda_{\rm obs}}{D_{\rm dish}}\right)^2 \left(\frac{1}{2 N_{\rm
dish} t_{\rm pix} \Delta \nu}\right) W_{\rm beam}^2(\ell)
\label{noise}
\end{equation}
In the above expression, $N_{\rm dish}$ denotes the number of interferometer
dishes, each assumed to have the diameter $D_{\rm dish}$, and $\lambda_{\rm
obs}$ is the observed wavelength. The  $\bar {T}(z)$ is the mean
brightness temperature at redshift $z$ defined by:
\begin{equation}
	\bar T(z) \simeq {44} \ {\mu {\rm K}} \left(\frac{\Omega_{\rm
HI}(z)h}{2.45\times10^{-4}} \right)\frac{(1+z)^2}{E(z)}
	\label{eq:tbar}
\end{equation}
where $E(z)=H(z)/H_0$ is the normalized Hubble parameter at that redshift.
The $T_{\rm sys}$ is the system temperature, calculated following $T_{\rm
sys} = T_{\rm inst} \ +  \ 60 \ {\rm K} \left(\nu/300  \ {\rm MHz}
\right)^{-2.5}$ where $T_{\rm inst}$ is the instrument temperature and $\nu$ is
the observing frequency.
The quantity $W_{\rm beam}^2(\ell)$ denotes the beam window function due to the finite angular resolution of the instrument operating in single-dish mode \citep[e.g.,][]{pourtsidou2016, pourtsidou2016a}, and is given by:
\begin{equation}
 W_{\rm beam}^2(\ell) = \exp\left[\frac{\ell (\ell + 1) \theta_B^2}{8 \ln 2}\right]
\end{equation}
where $\theta_B \approx \lambda_{\rm obs}/D_{\rm dish}$ is the beam FWHM of a single dish.

The integration time per beam is $t_{\rm pix}$ and the $\Delta \nu$ denotes
the frequency band channel width, which is connected to the tomographic
redshift
bin separation $\Delta z$. For the purposes of the noise calculation, we assume that
$\Omega_{\rm HI}(z)h = 2.45\times10^{-4}$, independent of redshift.

For the CHIME-like experiment we consider, we use the full interferometric noise expression, which is given by:
\begin{equation}
 N_{\ell, {\rm CHIME}}  = \frac{4 \pi f _{\rm sky}}{{\rm FoV}n_{\rm base}(u)n_{\rm pol} N_{\rm beam}  t_{\rm tot} \Delta \nu} \left(\frac{T_{\rm sys}}{\bar{T}}\right)^2
\left(\frac{\lambda_{\rm obs}^2}{A_{\rm eff}}\right)^2.
\label{noisechime}
\end{equation}

In the above equation, $N_{\rm beam}$ is the number of independent beams, and for CHIME, $N_{\rm beam} = N_f \times N_{\rm cyl}$ where $N_f = 256$ is the number of feeds, and $N_{\rm cyl} = 4$ is the number of cylinders.  Each is assumed to have the effective area $A_{\rm eff} = \eta L_{\rm cyl} W_{\rm cyl}/N_f$, where $\eta = 0.7$ , $W_{\rm cyl} = 20\,\mathrm{m}$ is the width of each cylinder, and $L_{\rm cyl} = 100\,\mathrm{m}$ is its length.  The total integration time, denoted by $t_{\rm tot}$ is taken to be 1 year for the CHIME-like survey considered here. The $n_{\rm pol}$ is the number of polarisation channels (taken to be 2). The baseline number density is $n_{\rm base}(u)$, expressible in terms of the multipole $\ell$ via $u=\ell/(2\pi)$. This quantity is approximated as independent of $u$ up to a maximum baseline length $u_{\rm max}$, viz. $n_{\rm base}(u) = N_{\rm beam}^2/(2 \pi u_{\rm max}^2)$. The $u_{\rm max}$ denotes the longest baseline $d_{\rm max}$ measured in wavelength units, $u_{\rm max} = d_{\rm max}/\lambda_{\rm obs}$,  with $d_{\rm max} = 269\,\mathrm{m}$ for CHIME \citep{obuljen2018}. The field of view for the CHIME interferometer is approximated as ${\rm FoV} \approx \pi/2 \times \lambda/W_{\rm cyl}$ \citep{newburgh2014} for this configuration.

}

The noise for the galaxy survey is taken to be the (Poisson) shot noise,
calculated as $N_{\ell, {\rm gal}} = n_{\rm gal, bin}^{-1}(z)$ where $n_{\rm
gal, bin}(z)$ is
the number density of galaxies per steradian in the bin centred at redshift
$z$.  Given the selection function of the galaxies, $\phi(z)$ defined in
\eq{selfunc}, this quantity is computed as:
\begin{equation}
 n_{\rm gal, bin} (z) = \int_{z - \Delta z /2}^{z + \Delta z /2} \phi(z') dz'
\end{equation}

Finally, the variance of the forecasted angular power spectrum is calculated as:
\begin{eqnarray}
&& (\Delta C_{\ell, \times})^2 = \frac{1}{(2 \ell + 1)  \Delta \ell f_{\rm sky,
\times}}  \nonumber \\
 &&\left[\left(C_{\ell, \rm HI} + N_{\ell, \rm HI} \right) \left(C_{\ell, \rm
gal} + N_{\ell, \rm gal} \right)
+  C_{\ell, \times}^2  \right]
\label{clerror}
\end{eqnarray}

In the above expression, the quantity $f_{\rm sky, \times}$ denotes the sky
coverage of the overlap between the surveys. For simplicity, an optimistic
complete overlap is assumed, and hence throughout this work, $f_{\rm sky,
\times}$ denotes the smaller of the two sky coverages of the galaxy and HI
redshift survey respectively. We use 15 $\ell$-bins between $\ell = 1$ and $\ell
= 1000$, logarithmically spaced with $\Delta \log_{10} \ell = 0.2$.

\begin{table}
\begin{center}
    \begin{tabular}{ | c | c | c | c | c | c | c | p{5cm} |}
    \hline
 Astrophysical &  & Cosmological & & \\
 \hline
 log ($v_{\rm c,0}$/km s$^{-1}$) & 1.56  & $h$ & 0.71  \\
 $\beta$ & -0.58 & $\Omega_m$ & 0.28 \\
 $Q$ &  4.6 & $\Omega_b$ & 0.0462 \\
 &  &  $\sigma_8$ & 0.81\\
  &  & $n_s$ & 0.963 \\
 \hline \\
    \end{tabular}
\end{center}
\caption{Fiducial values of the astrophysical and cosmological parameters
considered. {{Astrophysical parameters come from the best-fitting values of the
halo model for neutral hydrogen (Paper II), and that of the the galaxy
$Q$ parameter from the `blue5' galaxy sample in \citet{cresswell2009}. The
cosmological parameters are in good agreement with most available observations,
including the latest Planck results \citep{planck}.}}}
 \label{table:astrocosmo}
\end{table}

We use a Fisher forecasting formalism to place constraints on the cosmological
and astrophysical parameters, given the experimental configuration under
consideration. The Fisher matrix is computed as follows:
\begin{equation}
F_{ij} = \sum_{\ell} \frac{1}{(\Delta {C_{\ell, \times}})^2}\frac{\partial
C_{\ell, \times}}{\partial p_i} \frac{\partial C_{\ell, \times}}{\partial p_j}
\label{fisher}
\end{equation}
where the sum is over the range of $\ell$'s probed, and the $p_i$'s denote the
individual parameters.

The following parameters are used for the computation of the cross-power
spectrum of HI and galaxy surveys:
\begin{enumerate}
\item  The HI-based astrophysical parameters include $v_{c,0}$, $\alpha$, and
$\beta$ used in the $M_{\rm HI} (M)$ relation, and the two parameters $c_{\rm
HI, 0}$ and $\gamma$ for the HI profile,

\item the galaxy astrophysics contains the three parameters $b_0$, $b_{\rm
ls}$, $Q$ and $A$ used in the large-scale and the scale-dependent part of the
bias respectively, and

\item the cosmological parameters are the Hubble parameter $h$, the baryon
density $\Omega_b$, the spectral index $n_s$, the power spectrum normalization
parameter $\sigma_8$, and the cosmological matter density, $\Omega_m$.
\end{enumerate}

Of the HI astrophysical parameters, only two, viz.  the cutoff and the slope of
the HI-halo mass relation, i.e. $v_{c,0}$ and $\beta$ are relevant for
forecasting with HI intensity mapping surveys (see Paper I for details). While
we use all the galaxy parameters to model the bias for various surveys, we vary
only the parameter $Q$ encoding the scale-dependence of the bias.
Throughout the analysis, the cosmology adopted is flat, i.e. $\Omega_{\Lambda}
= 1 - \Omega_m$. The fiducial values of the cosmological and astrophysical
parameters are listed in Table \ref{table:astrocosmo}.

For calculating the standard deviations of the various cosmological and
astrophysical parameters, we use a procedure similar to Paper I: we consider
equal sized redshift bins of width $\Delta z \approx 0.5$ each, spanning the
desired
cross-correlation range in redshift, and evaluate the Fisher matrices $F_{ij}$
given by \eq{fisher} at the midpoints of each of the bins. The cumulative
Fisher
matrix for the $z-$range is derived from tomographic addition of the bins:
$F_{ij, \rm cumul} = \sum_{\Delta z \in z} F_{ij}$, which is the sum of the individual Fisher matrices, $F_{ij}$ in each of the $z-$bins of width $\Delta z$ contained between 0 and $z$.
From the cumulative Fisher matrix, the standard
errors in the parameters are computed for various cases. We ignore the
effects of cross-correlations between individual bins and those between galaxies
and HI in adjacent bins.

\section{Experiment combinations}
\label{sec:surveys}

For each galaxy survey, the specifications include the large-scale galaxy bias,
$b_{\rm gal, ls}$, the selection function $\phi(z)$ and the total number density
of galaxies, $n_{\rm gal}$. These as
well as the survey properties of the HI surveys are listed together in Table
\ref{table:surveyprops}.

 \begin{table*}
\begin{center}
    \begin{tabular}{ | c | c | c | c | c | c | c | c | c| p{1cm} |}
    \hline
   &   Galaxy & &  & HI &  &   Cross-correlation & & \\
    Configuration & $n_g$ {(\rm amin$^{-2}$)} & $b_{\rm ls, gal}$ & $T_{\rm
inst}$ (K) & $N_{\rm dish} $ & $D_{\rm dish} $ (m.)  & $f_{\rm sky, \times}$  &
$z_{\rm
bins}$\\ \hline
   CHIME-DESI  & 0.33  & 0.84/$D(z)$ & 50 & 1280 & 20  & 0.44 & [0.8, 1.2, 1.6]
\\
     MeerKAT-DES  & 8 & $(1.07 - 0.35 z)^{-1}$ &  29 & 64 & 13.5  &  0.12 &
[0.5, 1.0, 1.4] \\
     SKA I MID-LSST &  26 & $1.46(1 + 0.84 z)$ &  28 & 190 & 15  & 0.48 &
[0.082, 0.58, 1., 1.5, 2., 2.3, 3.06] \\ \hline
    \end{tabular}
\end{center}
\caption{Various experimental configurations considered in this work.}
 \label{table:surveyprops}
\end{table*}

\subsection{CHIME and DESI}

\begin{figure}
\begin{center}
\includegraphics[scale = 0.6, width = \columnwidth]{./selfunc_desi.pdf}
\caption{Redshift selection function constructed from the forecasted number
counts for DESI ELG galaxies in \citet{desi2016}, along with its fitted
functional form represented by \eq{desisel}.}
\label{fig:selfunc}
\end{center}
\end{figure}

The redshift coverage of the cross-correlation is $0.8 < z < 1.8$. The CHIME
autocorrelation survey runs over $z \sim 0.8 - 2.5$.
The DESI sample is assumed to correspond to the Emission Line Galaxy  (ELG)
survey \footnote{This is the largest sample of galaxies for which DESI will
obtain spectroscopic redshifts over $z \sim 0.6-1.8$.}, with the bias factor
$b_{\rm gal, ls} = 0.84/D(z)$ where $D(z)$ is the
growth factor.   The selection
function for DESI is constructed by numerically fitting to the number counts in
the ELG forecasts over $z \sim 0.6 - 1.8$, Table 2 of \citet{desi2016}. It is
found that this selection function can be modelled as:
\begin{equation}
 \phi(z) \propto (z/z*)^2 \exp(- (z/z_0)^{\beta})\
 \label{desisel}
\end{equation}
where $z* = 1.96, z_0 = 1.14$ and $\beta = 4.36$. This selection function, as
well as the raw number counts forecasted for DESI ELG galaxies, are plotted in
Fig. \ref{fig:selfunc}. The surface number density of
galaxies is $n_{\rm gal} \approx 0.33$ arcmin$^{-2}$ (corresponding
to roughly 1200 galaxies per square degree), which is consistent with the
estimates for the numbers of ELG targets in \citet{desi2016}.

\begin{figure}
\begin{center}
\includegraphics[scale = 0.6, width = \columnwidth]{./clexample.pdf}
\caption{Cross-correlation power spectra for a CHIME-DESI like survey, at the
three
redshifts of interest. Error bars indicate the expected standard deviation of
the angular power spectrum (\eq{clerror}) at
the lowest redshift ($z \sim 0.8$).}
\label{fig:clexample}
\end{center}
\end{figure}

The observing time $t_{\rm
pix}$ is assumed to be 1 hour (per pixel). The sky coverages for the individual
surveys are taken as $f_{\rm sky}$ (CHIME) =  0.61 (corresponding to $25000$
deg$^2$)
and $f_{\rm sky}$ (DESI) =  0.44 (corresponding to $18000$ deg$^2$) and the
DESI value is assumed for $f_{\rm sky, \times}$. We consider equal-sized
redshift bins of width $\Delta z = 0.5$ each.

Plotted in Fig. \ref{fig:clexample} are the cross-correlation angular power
spectra (computed following \eq{clcross}) for the
CHIME-DESI like configuration at three mean redshifts 0.8, 1.2 and 1.6. At the
lowest redshift, the error bars indicating the expected standard deviation on
the power spectrum (from \eq{clerror}) are also
plotted.

\begin{figure*}
\includegraphics[scale = 0.6, width =
\textwidth]{./barchart_evolution_chime_desi_nomarg.pdf}
\caption{Cross-correlation forecasts on astrophysics and cosmology for a
CHIME-DESI like survey. Fractional errors, $\sigma_A/A$ are plotted with $A =
\{h, \Omega_m, n_s, \Omega_b, \sigma_8, v_{c,0}, \beta, Q\}$, using information
from all the redshift bins available over the combined dataset. The left panel
shows the constraints on the cosmological parameters (i) without marginalizing
over astrophysics and (ii) marginalizing over astrophysics but including a
prior
based on the current knowledge of the astrophysical parameters. The right panel
plots the constraints on the astrophysical parameters for the case of marginalizing over cosmology but adding
the astrophysical prior. The extent of the astrophysical prior assumed is
plotted as the cyan band for each case in the right panel.}
\label{fig:forecasts}
\end{figure*}

In Fig. \ref{fig:forecasts} are plotted the cumulative fractional errors
(combining all the redshifts under construction) on the
forecasted cosmological and astrophysical (both HI and galaxy) parameters in
the following cases: (a) with fixed cosmology, i.e. without marginalization
over
the cosmological parameters, (b) with fixed astrophysics, and (c) marginalizing
over the galaxy and HI astrophysics, assuming a prior on the astrophysical
parameters coming from the current knowledge of the HI and galaxy data.  The
extent of these astrophysical priors are plotted in violet.  The HI
parameters are assumed to have the best-fit standard deviation values
constrained by the presently available data (see Table 3 of
\citet{hparaa2016}).{The galaxy
parameter $Q$ is taken to have a standard deviation of $1.78$, following the
discussion for
the `blue5' galaxy sample in \citet{cresswell2009}.}

The left panel  of Fig. \ref{fig:forecasts} shows that the constraints in the
CHIME-DESI cross-correlation case improve on the
 corresponding auto-correlation constraints using a CHIME-like
configuration alone, by factors of about 1.1-2 depending on the cosmological
parameter under consideration (comparing to Fig. 7 of Paper I). It is notable
that this improvement occurs even though the redshift coverage of the
cross-correlation is only about half that of the autocorrelation survey, and
illustrates the extent to which adding the galaxy survey information helps
improve
the cosmological constraints.

For all the three astrophysical parameters, there is a marked improvement on the current knowledge of the
astrophysics from the cross-correlation information, as represented by the
relative magnitudes of the violet and cyan/green bars. The  constraints improve by
factors
of 3.5, 3.1 and 3.3 respectively compared to their current priors.

\subsection{MeerKAT and DES}

\begin{figure*}
\includegraphics[scale = 0.6, width =
\textwidth]{./barchart_evolution_meerkat_des_nomarg.pdf}
\caption{Same as Fig. \ref{fig:forecasts}, for  a MeerKAT-DES like survey.}
\label{fig:meerkatdes}
\end{figure*}

\begin{figure*}
\includegraphics[scale = 0.6, width =
\textwidth]{./barchart_evolution_lsst_ska_nomarg.pdf}
\caption{Same as Fig. \ref{fig:forecasts}, for  an LSST-SKA like survey.}
\label{fig:lsstska}
\end{figure*}

 \citet{pourtsidou2016} discussed the potential for forecasting  lensing
convergence
parameters with a MeerKAT-DES survey in the southern hemisphere.\footnote{The
MeerKLASS
\citep[MeerKAT Large Area Synoptic Survey;][]{santos2017} proposes to
investigate galaxy evolution and cosmology using a 4000 deg$^2$ overlap with
the Dark Energy Survey (DES).} Here, we explore how such a cross-correlation
survey could potentially constrain the cosmological and astrophysical (both
galaxy bias and HI) parameters as in the previous case.

The redshift coverage of the survey is taken to broadly cover $z = 0.2$ to $z =
1.4$.
(Both DES and MeerKAT cover a similar redshift range, so the redshift overlap
is stronger between these two surveys.) The sky coverage is assumed to be all
of the DES survey, $5000 \  {\rm deg}^2$, thus assuming complete overlap. The
galaxy bias for the DES galaxies is given by the fitting formula
\citep{chang2016, pujol2016}: $b_{\rm gal, ls}^{-1} = 1.07 - 0.35z$.
The redshift selection function for the DES galaxies is taken to have the form
\citep[e.g.,][]{crocce2011}:
\begin{equation}
 \phi(z) \propto \left(\frac{z}{0.5}\right)^2
\exp\left(-\left(\frac{z}{0.5}\right)^{1.5}\right)
\end{equation}
The
surface number density of the DES galaxies under consideration is $n_{\rm gal}
=
$ 8 arcmin$^{-2}$, consistent with the estimates in \citet{becker2016}.

The fractional errors on the parameters considered are plotted in Fig.
\ref{fig:meerkatdes}. Errors in the 'fixed astrophysics case' remain essentially
unchanged from those with the astrophysical prior.
The
constraints on the cosmological parameters are of the same order of magnitude as those from the MeerKAT autocorrelation survey.
The astrophysical constraints improve over the current knowledge
of these parameters as
seen by the relative magnitudes of the cyan and the green bars, by
factors of about 1.5-2 for $v_{\rm c,0}, \beta$ and $Q$ respectively.

\subsection{SKA I MID and LSST}

The Large Synoptic Survey Telescope (LSST) survey parameters are taken to be
\citep{abell2009}, see also \citet{chang2013,ferraro2018}: (i) the galaxy
surface number density $n = 26$ arcmin$^{-2}$, and (ii) large scale bias
$b_{\rm gal, ls} = 1.46(1+ 0.84z)$. The redshift coverage is from $z \sim 0-3$
which spans both the SKAI-MID (B1 and B2) bands. The redshift selection function
of the survey is taken to be \citep{chang2013}:
\begin{equation}
\phi(z) \propto z^{1.28} \exp\left(-\frac{z}{0.41}\right)^{0.97}
\end{equation}
The sky coverage of LSST is
assumed to be 20000 deg$^2$ and that of the SKA-I is 25000 deg$^2$, and hence
the LSST coverage is used for calculating $f_{\rm sky, \times}$ (assuming
complete overlap).

This configuration leads to the tightest constraints on all the cosmological
and astrophysical parameters as shown in Fig. \ref{fig:lsstska}, with all
relative errors being about a few percent. It also leads to an substantial improvement in the
astrophysical constraints as compared with those from the present data.

The relative errors on the cosmological parameters reach values down to  $\sim
0.01$ with
this configuration. Constraints on $h$, $\sigma_8$ and $n_s$ with the
astrophysical prior improve by factors
of a few to ten, compared to the corresponding values from the SKA I -
MID like
autocorrelation survey alone (shown in Paper I), while those on
$\Omega_m$ and $\Omega_b$ improve by factors 2-5. Further, the astrophysical
parameters
improve over their current priors by factors of 2.7, 2.8 and 3.3 for $v_{c,0},
\beta$ and $Q$ respectively.

\section{Conclusions}
\label{sec:conclusions}
In this paper, we have explored combining upcoming HI intensity mapping surveys
with wide field galaxy optical surveys to improve available constraints on
astrophysical and cosmological parameters over $z \sim 0-3$ in the
post-reionization universe. Using the $\Lambda$CDM cosmological
parametrization,
a halo model framework for HI driven by currently available data, and available
optical
galaxy parametrizations, we have studied the extent to which these
constraints improve over their current uncertainties due to cross-correlation
measurements. We also note the improvement in the constraints compared to those
from the corresponding HI auto-correlation surveys alone.

For all three survey cases considered (a CHIME-DESI-like survey in the northern
hemisphere, and (ii) a MeerKAT-DES-like and LSST-SKA-like survey in the
southern
hemisphere), we find that the cross-correlation leads to
improvements in measurement of both astrophysical and cosmological parameters,
though the extent of improvement depends on the parameter under consideration. The significant benefit of cross-correlation (particularly in the MeerKAT-DES-like and CHIME-DESI-like configurations) lies in the improvement of astrophysical constraints.
The halo model framework allows us
to place realistic priors on the HI astrophysics from the currently available
data.
With the LSST-SKA combination, all the parameter constraints (both
astrophysical
and cosmological) reach levels below about 20 percent, even
without the assumption of cosmological priors.

The astrophysical forecasts for the HI and galaxy parameters improve
substantially (by factors of a few) over their current priors, with
the help of the cross-correlation measurements. This holds even in the presence
of the additional galaxy
parameter $Q$, which is seen to have comparable constraints though its prior
knowledge is assumed to be more uncertain.
The cosmological forecasts improve by factors of about a few (depending on the configuration) over those from
the corresponding autocorrelation surveys alone.

We note that the foregrounds, which may be the limiting systematic, are
excluded from the noise calculation in the present study. However, in the case
of cross-correlation measurements, the foregrounds for the two individual probes
are expected to be significantly
uncorrelated and thus lead to negligible effects \citep[as shown for the case
of
non-smooth foregrounds in 21 cm cross correlations with LBG surveys in,
e.g][]{navarro2015}. Recent studies \citep{foreground12019, foreground22019, breysse2019}
describe ways in which the signal, in the presence of foregrounds, may be
reconstructed from 21 cm intensity mapping data cross-correlated with other
tracers, such as the CMB or optical (both spectroscopic and photometric) galaxy
surveys. { A couple of caveats in this respect, however, are worth mentioning.

(i) First, we know that 21 cm experiments in autocorrelation suffer from bright foregrounds, which must be removed effectively, leading to a loss of low-$k$ modes (this effect is particularly prominent for the low-$k$ modes in the radial direction, which correspond to the largest scales). Hence, these Fourier modes are likely to be lost from the survey information, even if systematic biases caused by foregrounds can be disregarded.

(ii) We have not explicitly modelled the effects of the galaxy photo-$z$ errors in the present analysis. Since fairly broad redshift bins are used (with $\Delta z \sim 0.5$), the effect of these errors may be largely mitigated for the case of the galaxy-galaxy autocorrelation. However, if foreground filtering results in the removal of radial modes in the 21 cm surveys up to some maximum (as described in the previous point), then there may be little overlap between the 21 cm radial modes and the galaxy radial modes that remain for the cross-correlation. This effect may need careful treatment when foreground cleaning or avoidance is being considered.
}

Just as in the previous study with auto-correlation data alone (Paper II),  we
see
that the overlap in redshift coverage is extremely
important in tightening forecasts (due to more information coming from the
addition of independent tomographic bins). Also, choosing similar sky area
overlap between the HI
and galaxy surveys (presently assumed to have complete overlap)  would lead to
better constraints on the parameters.
Extending these approaches towards intensity mapping with other emission lines
(CO, CII) would enable
us to potentially form a comprehensive picture of galaxy evolution at the
scales
of the
ISM. Ultimately, combining both auto-and cross-correlation forecasts, possibly
with cosmological priors from present and future CMB experiments, would provide
the tightest possible constraints exploiting the synergy of CMB, HI and galaxy
surveys. This would be a powerful tool to explore more parameters in
cosmological models such as e.g., testing modifications to general relativity
at
the largest scales \citep[e.g.,][]{hall2013}  with future wide-field surveys.

\section{Acknowledgements}
We thank Stefano Camera, Ue-Li Pen and Ren\'{e}e Hlo{\v z}ek for useful
discussions, {{and the anonymous referee for a detailed and helpful report that improved the content and quality of the presentation.}} HP's research was supported by the Tomalla Foundation. AA is supported by the Royal Society Wolfson Fellowship.

\bibliographystyle{mnras}