Iterative map-making with two-level preconditioning for polarized Cosmic Microwave Background data sets
A worked example for ground-based experiments
Iterative map-making with two-level preconditioning for polarized CMB data sets
G. Puglisi et al.
Giuseppe Puglisi1<EMAIL_ADDRESS>Davide Poletti1   Giulio Fabbian  2,1 Carlo Baccigalupi 1,3 Luca Heltai1Radek Stompor4

SISSA- International School for Advanced Studies, Via Bonomea 265, 34136 Trieste, Italy
Institut d'Astrophysique Spatiale, CNRS (UMR 8617), Univ. Paris-Sud, Université Paris-Saclay, bât. 121, 91405 Orsay, France
INFN-National Institute for Nuclear Physics, Via Valerio 2, I-34127 Trieste, Italy
AstroParticule et Cosmologie, Univ Paris Diderot, CNRS/IN2P3,CEA/Irfu, Obs de Paris, Sorbonne Paris Cité, France

An estimation of the sky signal from  streams of Time Ordered Data (TOD) acquired by the Cosmic Microwave Background () experiments is one of the most important steps in the context of data analysis referred to as the map-making problem. The continuously growing data sets render the map-making problem progressively more challenging in terms of computational cost and memory in particular in the context of ground based experiments with their operational limitations as well as the presence of contaminants.
We study a recently proposed, novel class of the Preconditioned Conjugate Gradient (PCG) solvers
which invoke two-level preconditioners in the context of the ground-based CMB experiments. We compare them against the PCG solvers commonly used in the map-making context considering their precision and time-to-solution.
We compare these new methods on realistic, simulated data sets reflecting the characteristics of current and forthcoming ground-based experiment.
We develop an embarrassingly parallel, divide-and-conquer implementation of the approach where each processor performs a sequential map-making for a subset of the TOD.
We find that considering the map level residuals the new class of solvers permits achieving  tolerance of up to 3 orders of magnitude better than the standard approach, where the residual level often saturates before convergence is reached. This often corresponds to an important improvement in the precision of the recovered power spectra in particular on the largest angular scales. The new method also typically requires fewer iterations to reach a required precision and thus shorter runtimes for a single map-making solution. However, the construction of an appropriate two-level preconditioner can be as costly as a single standard map-making run. Nevertheless, if the same problem needs to be solved multiple times, e.g., as in Monte Carlo simulations, this cost has to be incurred only once, and the method should be competitive not only as far as its precision but also its performance is concerned.

§ INTRODUCTION

Over the last decades, several experiments have looked into the CMB polarization anisotropies aiming at  discovering a stochastic background of gravitational waves produced during the inflationary phase of our Universe encoded in the $B$-modes,  the divergence-free pattern in polarization. Indeed, the amplitude of the CMB $B$-mode polarization anisotropies at the scales larger than $1$ degree, conventionally parameterized with a tensor-to-scalar ratio, $r$, is thought to be directly related to the energy scale of inflation ($\sim 10^{16}$ GeV). These primordial $B$-modes have not been detected yet and further progress in both the control of the diffuse polarized emission from our Galaxy (involving widely the microwave frequency regime [Ade et al., 2015]) and in the sensitivity of the experimental set-ups is necessary in order to reach such a goal. At the sub-degree angular scales,  $B$-modes are produced by the gravitational lensing due to  large scale structures intervening along the photon path travelling towards us.  Evidence for these lensing $B$-modes was first provided via cross-correlation of the CMB polarization maps with the cosmic infrared data [Hanson et al., 2013, The Polarbear Collaboration et al., 2014] and via constraining the small-scale $B$-mode power [The Polarbear Collaboration et al., 2014] and they have been since then characterised with increasing accuracy [Louis et al., 2017, Keisler et al., 2015, BICEP2 Collaboration et al., 2016, The Polarbear Collaboration et al., 2017].
While these past experiments have observed the microwave sky with arrays of thousands of detectors often focusing on small sky patches, the forthcoming experiments are planned to observe bigger patches with at least tens of thousands of detectors, producing as a result, Time Ordered Data (TOD) including tens and hundreds of billions of samples.

The recovery of the sky signal from these huge, noisy time streams, a process called map-making, represents one of the most important steps in data analysis and, if the detector noise properties and scanning strategy are known,  map-making becomes a linear inverse problem.
The Generalized Least-Squares (GLS) equation provides an unbiased solution to map-making for an arbitrary choice of weights given by a symmetric and positive definite matrix [Tegmark, 1997].
Moreover, if we consider as the weights the inverse covariance of the time-domain noise, the GLS  estimate is also a minimum variance and a maximum likelihood solution to the problem. However, a computation of the solution in such a case may require either an explicit factorisation of a huge, dense matrix [Tegmark, 1997, Borrill, 1999, Stompor et al., 2002]  or an application of some iterative procedure [Wright, 1996, Oh et al., 1999, Doré et al., 2001, de Gasperis et al., 2005, Cantalupo et al., 2010]. These latter involve typically several matrix-vector multiplications at each iteration step.
What makes the map-making problem particularly challenging are the sizes of the current and forthcoming CMB data sets which are directly related to the number of floating point operations (flops) needed to achieve the solution and to the memory requirements due to the sizes of the arrays required for it.
Both these factors set the requirements on computational resources and indeed many current data analysis pipelines opt for  massively parallel computing platforms. However, even in such circumstances efficient algorithms are necessary to ensure that the analysis can be indeed performed. The computational complexity of the algorithms involving  an explicit matrix inversion is $ \mathcal{O}(N_p ^3)$ flops,  where $N_p$ is the  number of pixels in the map, and therefore they are suitable only for the cases when the estimated sky maps do not involve many sky pixels.
However, whenever feasible the direct approaches can yield high-precision, unbiased estimates of the sky signal  <cit.>. However, the next generations of the ground experiments, -Stage III  [Arnold et al., 2014, Henderson et al., 2016, Benson et al., 2014] and IV [Abazajian et al., 2016]), are expected to observe significant fractions of the entire sky with high resolution and thus resulting in maps with  $N_p \simeq {\cal O}(10^6)$,  rendering the direct approaches prohibitive even for the largest forthcoming supercomputers.

In this context iterative methods have offered an interesting alternative.

They involve algorithms within the class of  Krylov methods <cit.>, which avoid the explicit inversion of the linear system matrix by constructing an approximate solution which is iteratively improved on. The  computational complexity of such methods is mostly driven by  matrix-vector products, which need to be performed repeatedly on each iteration. These require at most ${\cal O}(N_p^2)$ flops and can be performed at much lower cost in the specific case of the CMB map-making (see Sect. <ref>), where such matrix-vector products can be computed matrix-free, i.e., without ever assembling explicitly the system matrix in  memory [Cantalupo et al., 2010]. To date, most of iterative solvers have been based on the  Conjugate Gradient (CG) method applied to a preconditioned system of map-making equation and involved a simple block-diagonal preconditioner (see Eq. (<ref>)). While these solvers performed usually very well <cit.>, the anticipated data sets motivate a search for better, more efficient algorithms [Grigori et al., 2012, Naess et al., 2014, Szydlarski et al., 2014, Huffenberger & Næss, 2018].

In this paper we apply the methodology proposed by Szydlarski et al., 2014 to a reconstruction of maps from simulated data of a modern, ground-based experiment. This new class of approaches involves constructing and applying a more involved, two-level preconditioner. Our simulations are informed by the experiences derived from the deployment and analysis of the experiment whose results from the first   two seasons of data have been recently published in The Polarbear Collaboration et al., 2017. In Sect. <ref>, we briefly introduce the formalism of the  map-making problem in presence of time domain filtering operators typical for a ground based experiment. In Sect. <ref>, we describe the iterative approach and the two different methodologies adopted in the analysis. In Sect. <ref> we further describe how the filters introduce degeneracies on the estimation of the maps. The main results of this paper are presented in Sect. <ref>, where we show a close comparison of the performances of the commonly-used block-diagonal and the two-level preconditioners applied onto realistic simulated data sets. These include signal as well as sky coverage and noise properties representative of the forthcoming generation of polarization experiments. Finally, in Sect. <ref> we summarize our main results and conclusions.

§ MAP-MAKING IN GROUND-BASED EXPERIMENTS

The input data of the map-making procedure are the calibrated TODs collected in a single time-domain vector $d$ of size $N_t$ containing all measurements performed during a certain period of time by all the detectors of a CMB experiment. The measurements can be modelled as the sum of an astrophysical signal and  measurement noise, ${n}_t$. The astrophysical contribution to a measurement taken at time $t$ is given by the sky signal, $s_p$, in pixel $p$ observed at time $t$ and which is already convolved with the instrument response, assumed hereafter to be axially symmetric. The correspondence between the sky pixel, $p$, and the time, $t$, can be encoded by  a sparse and tall ($N_t \times N_p$) matrix, ${P}_{tp}$.  The data model can be then written as:

\begin{equation}
{d}_t= {P}_{tp}  {s}_p + {n}_t,
\label{eq:probl}
\end{equation}

or in the matrix form as,

\begin{equation}
d= Ps + n. \label{eq:datamod}
\end{equation}

Here, $s$ stands for the map to be estimated.

The structure of the pointing matrix encodes the scanning strategy of the experiment and depends on whether the detectors are sensitive or not to the polarization. In the former case, the sky signal is described by three Stokes parameters $I,\,Q,\,U$ in every pixel $p$ of the map, i.e., $s_p=(I_p,Q_p,U_p)$, and a measurement by a polarization-sensitive experiment taken at time $t$ can be written explicitly as,

\begin{equation}
d_t= I_{p_t} +Q_{p_t}\cos(2 \phi_t) + U_{p_t} \sin (2 \phi_t) + n_t, \label{eq:datamodel}
\end{equation}

where $\phi_t$ is the angle of the detector projected onto the sky coordinates at time $t$. The pointing matrix has in this case three non zero entries per row.
We further assume the noise to have vanishing mean $\langle n \rangle =0 $ and defined by the noise covariance matrix $\mathcal{N}$.

Under these assumptions the map-making is a linear inverse problem of estimating the sky signal, $s$, from the data, $d$, given the data model as in Eq. <ref>. This is a linear statistical problem whose solution is provided by a GLS,

\begin{equation}
\hat{s} = (P^\dag   W P)^{-1} P^\dag   W d,  \label{eq:estimator}
\end{equation}

yielding an unbiased  estimator [Tegmark, 1997] for any choice of a positive definite matrix $W$.  In particular, if $W= \mathcal{N}^{-1}$ and the noise is Gaussian distributed the estimator in Eq. (<ref>) becomes minimum variance.

§.§ The filtering operator

The raw TODs are often contaminated by some unwanted signals that are not astrophysical in their origin, such as the ground pickup or the atmospheric contributions, or their noise properties display strong, long-term correlations commonly referred as $1/f$ noise. All these contributions are usually filtered out from the data.

In such cases the template of the unwanted signal, $T$, is known while its amplitude $y$ is not. What a filtering operation is required to do is to remove a component of the TOD contained by the subspace spanned by the columns of $T$,

\begin{eqnarray}
d^{\prime}\ \ &\equiv& (\mathbb{1} - T(T^\dag   T)^{-1} T^\dag  ) d = F_T d \\
F_T T&=&0,\label{eq:filt}
\end{eqnarray}

so that $d' \cdot T_i  = 0$, for any template $T_i$ included as a column of the  template matrix, $T$. The most general form of the filtering operator involves also weighting by a full-rank weight matrix, $W$, and reads,

\begin{equation}
d'\equiv (W - WT (T^\dag   WT)^{-1} T^\dag  W) d=F_T d.
\label{eq:filtnoise}
\end{equation}

With the above definition of $F_T$,   it is therefore possible to generalize (<ref>) to  [Poletti et al., 2017]

\begin{equation}
\hat{s} = (P^\dag  F_T P)^{-1} P^\dag   F_T d. \label{eq:unb}
\end{equation}

Notice that the filtering operator does not change the properties of the estimator in (<ref>). It is still unbiased,

⟨ŝ -s ⟩= ⟨(P^†F_T P)^-1 P^†F_T n  ⟩=0,

and if we consider $W = \mathcal{N}^{-1}$,  it is minimum variance.

§ PRECONDITIONED ITERATIVE SOLVERS

We can rewrite (<ref>) as a linear system,

\begin{align}
\left( {P}^\dag    {F_T P}\right) \hat{s}&=  {P}^\dag   { F_T d}, \label{eq:linearsystem} \\
&\Downarrow \nonumber \\
\end{align}

where ${A}$ is a symmetric and positive definite (SPD) matrix.

The CG algorithm is particularly attractive for large sparse or structured systems since it references the system matrix  only through its multiplication of a vector.
The convergence rate of the CG depends on the condition number of the system matrix, $\kappa$, [Golub & Van Loan, 1996], defined as the ratio of the largest to the smallest eigenvalue of a matrix.

To reduce the condition number of ${A}$,  a preconditioner matrix $M_P$ is applied to the linear system so that the condition number of matrix $M_P A $ is smaller. If this is the case, the CG  converges within a smaller number of iterations. This new algorithm is commonly referred as Preconditioned  CG (PCG) as it solves the preconditioned linear system,

\begin{equation}
M_P {A} x= M_P b.
\label{eq:pcg}
\end{equation}

It can be shown  [Golub & Van Loan, 1996] that the PCG convergence rate is strictly related to the condition number of the preconditioned matrix $M_P A$. In fact, after $k$ iterations of the PCG, the magnitude of the error is

\begin{equation}
\parallel e^{(k)} \parallel  \equiv\frac{\parallel x -x^{(k)} \parallel}{\parallel x \parallel} \leq \kappa(M_P A) \frac{\parallel r ^{(k)}  \parallel}{\parallel b \parallel},
\label{eq:errnorm}
\end{equation}

where $x$ is the true solution to (<ref>) and

\begin{equation}
r ^{(k)}\equiv  b - M_P A x ^{(k)} \label{eq:residuals}
\end{equation}

is the PCG residual at the $k$-th step.

§.§ The Jacobi  Preconditioner

From (<ref>) we can define the Jacobi Preconditioner:

\begin{equation}
\mbd^{-1} \equiv P^\dag \text{diag}(\mathcal{N}^{-1})   P.
\label{eq:mbd}
\end{equation}

This is not only trivial to compute, store and apply to a vector, but it also accounts for some of the eigenstructure of the actual system matrix, which is due to the inhomogeneity of the sky observations. These properties justify why $M_{BD}$ is the most popular and very successful preconditioner used in the current CMB map-making practice. We call it either the block diagonal or Jacobi preconditioner[We note that a more typical definition of the Jacobi preconditioner, i.e., $\text{diag}(A)$, would not have the same attractive computational properties because its computation would require handling a dense time domain square matrix, $\mathcal{N}$. The upside of the block diagonal preconditioner is precisely that it takes care of the scanning strategy-induced increase in the condition number without dealing with the complexity of the time domain processing.].

The effect of Jacobi preconditioners onto the eigenspectrum of   is to shift the largest eigenvalues  towards the unity, thus potentially decreasing the condition number of the preconditioned system. However, the nearly singular eigenvalues due to the noise correlations or the filtering will not in general be accounted for.
These are common for ground based experiments and consequently the convergence of the PCG with the block-diagonal preconditioner is often found unsatisfactory. Indeed, in extreme, albeit not uncommon, cases this is manifested as a saturation of the residuals level and lack of the actual convergence down to a required threshold <cit.>.

§.§ Two-level Preconditioners

An alternative preconditioner may be found among the class of the so called Deflation preconditioners  that have proven to be successful in presence of few isolated extremal eigenvalues.
They act as de-projectors from the so called deflation subspace, $\mathcal{Z}$. This subspace is   generated by $r$ linearly independent eigenvectors  related to the smallest  eigenvalues and constitute the columns of the deflation  matrix $Z$. This matrix is needed to define the projector $R$

\begin{equation}
R = \mathbb{1} - {A}Z(Z^\dag   {A}Z)^{-1} Z^\dag.
\end{equation}

The projector $R$ is $A$-orthogonal to any vector $ w \in \mathcal{Z}$ since $R{A}Z=0$.
In the exact precision algebra, $R$ would be a very efficient preconditioner, as for each steps of an iterative CG-like solver would be orthogonal to the null space of the $R{A}$. However, we deal with finite precision arithmetic and the zero eigenvalues are often as bothersome as the small ones due to the numerical precision of the machine.

This issue can be solved by combining the operator $R$ with the Jacobi preconditioner  as it has been proposed in Szydlarski et al., 2014

\begin{align}
\md\equiv&  \mbd R + ZE^{-1} Z^\dag   \nonumber \\
=& \mbd(\mathbb{1} - {A}Z(Z^\dag   {A}Z)^{-1} Z^\dag  ) + ZE^{-1} Z^\dag  , \label{eq:m2l}
\end{align}

where $E$ is the coarse operator, defined as $E = Z^\dag   {A}Z$.
$\md$ is referred as the two-level preconditioner and we note that it indeed fixes the issue of the  zero eigenvalues since they are rescaled all to one. Indeed,

\begin{equation}
\md {A}Z =Z.
\label{eq:ident}
\end{equation}

The dimension of  the deflation subspace, given by $dim (\mathcal{Z})= r$ is by  construction much smaller than $N_p$, and it is straightforward to invert the matrix $E$.  Moreover, as  is SPD, so is $E$.

We can summarize the action of the $\md $ preconditioner,  applied on a vector $v$, as  a projection of the vector $v$ onto two  different  subspaces, namely $\mathcal{Z}$ and its orthogonal complement  $\mathcal{Y}$.  The components of $v$  are projected onto  $\mathcal{Z}$  via the $ZE^{-1} Z^\dag $ term in (<ref>). In this subspace, the inverse  is  very well approximated by $\md$,  since we have that (<ref>) holds for any $z \in \mathcal{Z}$. On the other hand,  $\md$ acts onto a generic vector $y \in \mathcal{Y}$ in the same way $\mbd$ does, since

Ay = Ay.

Thus, once $\md $ de-projects from the deflation subspace, it performs   the  PCG by means of the standard  preconditioner and it converges faster  since $\mbd {A}$  has a smaller condition number $\mathcal{O}(10)$ (due to the considerations made at the end of Sect. <ref>).

It may appear that in order to build the deflation subspace, one would require the knowledge of the entire eigenspectrum of  to determine the eigenvectors with the smallest eigenvalues.  However,  Szydlarski et al., 2014 has proposed that approximated eigenpairs derived with the help of the so called Ritz approximations (see Appendix <ref>) is sufficient for this purpose.

§ THE CASE OF THE GROUND BASED EXPERIMENTS.

A ground-based experiment, scanning the sky with a focal plane including thousands of polarization sensitive pixels, has  to cope with both atmospheric and ground emissions, which have to be treated on the time domain level. This can be achieved by applying filtering to the data as discussed in Sect. <ref>. The specific templates often applied in this context <cit.> are as follows.

§.§  Atmospheric emissions and noise correlations

Power Spectral Density of summing and differencing the signal from simulated data. Notice the $f_{knee}/f$ dependence at small frequencies, and the flattening due to white noise  above $1 \si{. \hertz}$. Solid blue (orange) line refers to a signal with a $f_{knee}= 1\, (0.05) \si{. \hertz}$.

Both the atmosphere fluctuations, $a$, as well as detector noise, $n$, introduce contributions correlated on long time scales. While such effects could be potentially suppressed by adopting an appropriate weight matrix, $W$, in practice such a solution is prohibitive given the sizes of the current and anticipated data sets in the time domain. In such cases the diagonal weight matrices, while straightforward to operate on, will  typically lead to poor quality estimates of the sky signal with strongly correlated, spurious features appearing along the scan directions.

Such long temporal modes can be however well approximated by an arbitrary linear combination of piece-wise low order polynomials. Collecting these templates in a matrix $B$, we can express the resulting residual as,

\begin{equation}
w \equiv  a + n - Bx.
\label{eq:polyfiltering}
\end{equation}

Since $Bx$ is supposed to capture the correlated component of $a+n$, $w$ is expected to be approximately described by a white noise statistics. Consequently, a diagonal weight matrix can be now sufficient for obtaining good quality sky maps. In our simulations we assume this approximation to be exact and thus simulate TODs with piecewise-stationary white noise.

Filtering these particular modes results  in removing from both noise and  signal long term trends present in the TODs, whose signal-to-noise ratio is usually very low. Even though, the stripes in the reconstructed map disappear we have to remember that the constraints on the large angular scales are weak. The system matrix $A$ encodes this information: the presence of the filtering operator (see the left-hand side of <ref>) results in low eigenvalues corresponding to long modes.

§.§ Ground pickup

Though ground-based experiment are designed to have very low far side lobes of the beam, the signal from the ground is not negligible compared to the one.  The elevation is typically constant during an observation period and therefore the ground signal can be considered as a function of  the azimuth. If we neglect contributions from other signals, the TOD data model can  be written as:

d= Ps + Gg +n.

Intuitively, we can think of the second term as  the ground template map $g$ projected to the time domain by means of the “ground-pointing matrix” $G$. This matrix has a column for each azimuthal bin, the entries of the column are equal to 1 whenever the azimuth of the pointing direction falls within the bin range, and they are zero elsewhere.

§.§ Map-making for the ground based experiments.

The effects discussed earlier in this Section have typically amplitudes significantly higher than those of the sky signals, which moreover do not average out efficiently while projected on the sky. For this reason they need to be treated explicitly in the map-making process by introducing appropriate filters.
As elaborated on in Poletti et al., 2017, the required filters, while dealing efficiently with the unwanted contributions, may however render the system matrix  ill conditioned, i.e. $\kappa\gg 1 $, implying the existence of degeneracies between a certain sky signal  $\tilde{s}$ and the amplitude of a certain template $\tilde{y}$. This means that some particular mode of sky signal is impossible to reconstruct  whenever the template $\tilde{y}$ is filtered out, since $P\tilde{s} = T \tilde{y} $ and $P^{\dag}F P\tilde{s} = P^{\dag} FT \tilde{y}=0$.

This may be particularly acute in the case of the ground pick-up filtering. For a constant-elevation one-detector observation filtering ground-stationary signal results in unconstrained modes that are constant in the right-ascension direction.  An extensive discussion on the properties and structure of the $F_T$ matrix can be found in Poletti et al., 2017.

Accumulating multiple detectors and observation periods can partially break these degeneracies, but the constraints on these modes will typically be weak.

As pointed out earlier, the presence of small-eigenvalues in the eigenspectrum of the system matrix, $A$, can significantly affect the convergence of the iterative solvers and can not be accounted for by the standard, block-diagonal preconditioner.

§ THE SIMULATED DATA SET

In this section we describe the  experimental setup we adopted to perform map-making runs.
We exploit the simulations capabilities of the  package[https://github.com/JulienPeloton/s4cmb] to produce simulated data sets for different experimental configurations of a ground-based experiment located in the Atacama desert in Chile at an altitude of  $\si{\num{5190}}~\si{\meter}$. We remind that the site location has implications for the properties of the observation. For example, since every pixel is observed at different elevations, the projection of the scan on the sky crosses the pixel with different direction, increasing the so-colled cross-linking (i.e. the coverage in the orientation of the attach angle).

We consider a 60-wide focal plane hosting dual-polarization pixels  sensitive to 148 with a fractional bandwidth of $26 \%$. The resolution of the telescope is assumed to be $3.5$ arcmin.

We consider three cases that differ for the target sky area and the sensitivity of the instrument, they are summarized in table <ref>.
The configurations labelled Small Patch (SP) and Large Patch (LP) refer to the characteristics of current and forthcoming experiments observing either small ($f_{sky} < 1\%$) or wide ($f_{sky} \gtrsim 1\%$) sky patches.  For both cases, the Noise Effective Temperature (NET) per detector is $\sim 500 \mu\mathrm{K \, \sqrt{s}}$  but  in LP we increase of about one order of magnitude the number of detectors in the focal plane, from $600 $ to $8,000$.

We  consider  an additional case to reproduce a next generation of ground-based experiments that will observe a wider fraction of the sky ($f_{sky} \sim 20 \%$) with an increased detector sensitivity NET$\sim 360 \mu\mathrm{K \, \sqrt{s}}$ and a larger number of detectors ($50,000 $). We refer to this setup as the Very Large Patch (VLP).

The simulated observations are divided into constant elevation scans (CESs) during which the telescope scans back and forth in azimuth at a speed of  $\si{\ang{0.4} \second^{-1}} $ and at constant elevation  (hereafter, we commonly refer to  each azimuthal sweep as  a subscan).   When the patch has moved out of the field of view, the telescope moves the elevation and azimuth toward the new coordinates of the patch and a new CES starts.
The samples are acquired at a rate  of 8, which given our scanning speed is sufficient to reach $\ell \sim 1200$. The number of samples per CES depends on the width of the subscan  and on the number of detectors performing the measurements.

Properties of the different scanning strategies

Case          $N_t$      $ N_p$    $f_{sky}$   $ NET_{array} \, [\mathrm{\mu K \sqrt{s}} ] $    Observation time [yr]

Small  Patch      $\si{\num{3d10}}$	   $\si{\num{4d4}}$    $0.1 \% $    20.4        2

Large  Patch      $\si{\num{3d10}}$      $\si{\num{2d6}}$    $5 \% $     5.6    1

Very  Large Patch     $ \si{\num{3d10}}$    $ \si{\num{1d7}}$    $ 20\%$     1.6    1

Using the simulated scanning strategies we scan an input map computed with the  routine of  the Hierarchical Equal Area Latitute Pixelization  (HEALPIX library, Górski et al., 2005) [http://healpix.sourceforge.net] and then add a white noise realization corresponding to the sensitivity of each experimental configuration. The input signal power spectrum has been computed with the CAMB package [Lewis et al., 2000] assuming the  2015 best fit cosmological parameters [Planck Collaboration et al., 2016] and $r=0.1$. Since we consider maps at the resolution of $3.5$ arcmin, we sampled the input CMB signal on a grid with $1.7$ arcmin resolution  (corresponding to a resolution parameter $\texttt{nside}=2048$).
The definition of the observed pixels is performed prior to the map-making procedure and it is based on discarding those samples which do not observed pixels with enough redundancy.

§.§ Data model

The simulated TODs acquired by each detectors can be expressed as in Eq. (<ref>). However, since the detectors are grouped in pairs sensitive to orthogonal polarization states of the light, commonly referred as $d_{top}$ and $d_{bottom}$, the signal coming from the two  detectors of a given pair can be combined in order to disentangle the total intensity and polarization signals without any loss of accuracy by summing and differencing two signal:

d^±≡1/2( d_top ±d_bottom).

Hence, one can independently estimate intensity and polarization (expressed via the $Q,U$ Stokes parameters) maps. Two separated data models can be written for the signal and noise component of the time streams:

\begin{align}
d^+_t \equiv&P_{tp}I_p +  n_t^+, \label{eq:d+}\\
d_t^- \equiv &P_{tp} ( \cos (2\phi_t) Q_p  + \sin(2\phi_t) U_p ) + n_t^-\label{eq:d-}
\end{align}

where $n^{\pm} $ is the noise term and can be analogously defined as $d^{\pm}$.   We note that in the following we will neglect possible systematic effects that can leave a residual intensity component in $d_t^-$, such as a calibration mismatch between the two detectors of a pair. The noise properties of the sum and difference time streams are different, as depicted in Fig.  <ref>, and therefore require a different set of low order polynomials to be filtered out. In our analysis, the time stream $d^{-}$ if filtered  by  zeroth and first order polynomials, $d^+$ by the first three order polynomials  and we assume that these filters completely remove the $1/f$ component. For this reason our noise simulation contain solely white noise. Simulating the correlated noise component is important when evaluating the end-to-end performances of an experiment, but in this paper we focus only on the performances of the map-making solver, which depend mostly on the scanning strategy and data processing adopted.
The sum-difference approach and the fact that $n_t^+$ and $n_t^-$ are uncorrelated allow to separate the intensity and the polarization reconstruction, we take advantage of this by focusing only on the latter for the rest of this paper.

The ground template is the same for summed and differenced data, and its column number is the same as the number of azimuthal bins ($100\div 1000$, depending on the width of the patch). Each azimuth bin has a fixed width of $\si{\ang{0.08}}$. The rows are as many as the number of samples in each CES $N^{CES}_t$.

For simplicity, in the following analysis, we do not build $F_{[G,B]}$, instead, we avoid the burden of explicit orthogonalization of the filters by using as the filter $F_T$ in Eq. (<ref>), a simplified filter given by $F_B W^{-1} F_G W^{-1} F_B$, which is explicitly symmetric and would have been equivalent to $F_{[G,B]}$, were the filters, $F_G$ and $F_B$, be orthonormal from the outset [Poletti et al., 2017].

§ CONSTRUCTING TWO-LEVEL PRECONDITONER

Ritz Eigenvalues $\mbd A$ estimated with different Arnoldi algorithm tolerance , $\epsilon_{Arn}$. An histogram version of this plot can be found in Fig. <ref>.

A construction of the two-level preconditioner requires knowledge of the deflation operator, ${\cal Z}$. We estimate it following the prescription of Szydlarski et al., 2014, which employs the Arnoldi algorithm to compute approximate eigenpairs of the matrix $B=\mbd {A}$. A suitable selection of these is then used to define the deflation operator, $Z$. This process has two free parameters, $\epsilon_{Arn}$ and  $dim(\mathcal{Z})$, that we discuss in the rest of this section and fix them in Sect. <ref> using numerical experiments.

The Arnoldi algorithm iteratively refines the approximate eigenpairs of the provided matrix and ends the computation when a given tolerance, $\epsilon_{Arn}$, is reached (see Appendix <ref> for more details). The lower the Arnoldi tolerance, the larger is the rank of the approximated $B$ and, consequently, the larger is the number and the accuracy of the estimated eigenpairs. In Fig. <ref>, the approximated eigenvalues are reported for several choices of $\epsilon_{Arn}$ and some specific choice of the system matrix corresponding to the small patch case discussed later. It shows that not only the number but also the range of the eigenvalues increase with smaller tolerance. This is intuitively expected since the Arnoldi algorithm relies on the power method <cit.>, and thus it first estimates the largest eigenvalues before moving to the smaller ones. Once the tolerance is as small as $\num{d-9}$ the range of the estimated eigenvalues starts to saturate. If we attempt to reach a threshold smaller than $\sim$d-12, the Arnoldi iteration proposes a new search direction that, due to the finite numerical precision, has no component linearly independent from the previous ones. Consequently, from that moment on the algorithm keeps producing eigenvalues equal to zero that are not eigenvalues of $B$ but just the sign that the Arnoldi algorithm has converged and exhausted its predictive power. In the studied cases we find that this typically requires  $\sim 150$ iterations.

The computational time required for each Arnoldi iteration is similar to the CG but the memory consumption can be very different: while the  Arnoldi needs space for as many vectors as the iteration number, the CG requires only few vectors in the memory regardless of the iteration number. However, this is not a problem if the size of the map is negligible compared to the time streams – a condition likely to be met in forthcoming ground-based experiments.

The other parameter in the construction of the preconditioner is the dimension of the deflation space. For any given Arnoldi tolerance, this can be either fixed directly by defining the number of the smallest eigenvalue and eigenvectors retained to construct ${\cal Z}$ or by defining a threshold below which the eigenvalues and eigenvectors are retained, $\epsilon_{\lambda}$.

(left)  PCG residuals  with $\md$ and $\mbd$ preconditioners for several choices of the size of the deflation subspace, Ritz eigenpairs are computed up to a fixed Arnoldi tolerance $\epsilon_{Arn}$. (right) The $\md$  is built by selecting  a fixed  number of Ritz  eigenvectors,   $dim (\mathcal{Z})=25$, computed by running the Arnoldi algorithm with several choices for $\epsilon_{Arn}$.

§ RESULTS AND DISCUSSION

In this section, we present the performance comparisons of the standard block-diagonal and the two-level preconditioners both applied on simulated noise or signal-only dataset observing with the scanning strategies listed in Table <ref>. Moreover, we focus on the reconstruction of the polarization component of the sky,  but the results for total intensity are similar and are  reported in Appendix <ref>.

§.§ Comparison methodology

We use three types of metric in order to estimate the level of accuracy achievable by each considered approach. First we consider the norm of the standard map-level residuals  as defined  in Eq. (<ref>) normalized by the norm of the right hand side,

\begin{equation}
\parallel \hat{r}^{(k)} \parallel \equiv \frac{\parallel r^{(k)} \parallel }{\parallel b\parallel }.
\end{equation}

This measure of convergence is naturally provided in the CG algorithm and, most important, does not require the knowledge of the true solution. It is indeed the one typically employed in real applications for measuring the reconstruction quality.

In order to get further insights, in this paper we also consider metrics that require knowledge of the exact solution, which is available only in the case of signal-only simulations.
We make use of the norm of difference between the true and recovered map,  defined in Eq. (<ref>), and the bin-by-bin difference between the power spectrum of the input map and reconstructed map, binned using equally spaced bins in multipoles, $\ell_b$,

\begin{equation}
\Delta \mathcal{C}_{\ell_b}^X \equiv \frac{|\mathcal{C}_{\ell_b}^{X,in} - \mathcal{C}_{\ell_b}^{X,out} |}{\sigma^X_{CV}} , \label{eq:deltapowsp}
\end{equation}

with $X= E, B $;  the differences are normalized with respect to the   cosmic variance  of the input map,

\begin{equation}
\sigma^X_{CV}(\ell_b,\mathcal{C}^X_{\ell_b} ) \equiv \sqrt{\frac{2}{(2 \ell_b +1) f_{sky } \Delta \ell_b}} \mathcal{C}^X_{\ell_b}.
\end{equation}

This power spectrum difference enables to check which scale in the maps are better constrained, and the normalization gives an estimate of how much the signal intrinsically  fluctuates. We stress that we compare against the power spectrum of the input map, not the power spectrum used to simulate it. Therefore, the normalization is just a reference value and $\Delta \mathcal{C}_{\ell_b}^X$ has no cosmic variance.

As the considered sky patches cover only a fraction of the sky, the power spectra are computed using a pure-pseudo power spectrum estimator  [Grain et al., 2009].  This is a pseudo power spectrum method [Hivon et al., 2002] which corrects the E-to-B-modes leakage  arising in presence of incomplete sky coverage [Smith & Zaldarriaga, 2007, Bunn et al., 2003, Lewis et al., 2001].

§.§ Setting the two-level preconditioner

We use numerical experiments to show the role of the two free parameters involved in the computation of the two-level preconditioner, $\epsilon_{Arn}$ and ${\rm dim}\,{\cal Z}$.
A sample of the results is shown in Fig. <ref>. In the left panel we fixed the Arnoldi tolerance to $\epsilon_{Arn}= 10^{-5}$ and change the dimension of the deflation space from ${\rm dim}\,{\cal Z} = 2$ to $50$, which corresponds to varying $\epsilon_\lambda$ from $0.01$ up to $1$.

The size of the deflation subspace affects strongly the steepness of the initial convergence. This is expected because, if we use all the vectors produced by the Arnoldi algorithm to construct the deflation subspace, the residual after the first iteration is related by construction to the Arnoldi tolerance. On the contrary, the case with the block-diagonal preconditioner corresponds to ${\rm dim}\,{\cal Z} = 0$. The more we include vectors in the deflation subspace, the more we approach ${\rm dim}\,{\cal Z} = 50$, which retains nearly all the vector produced by the Arnoldi iterations and indeed jumps immediately to a residual close to the Arnoldi tolerance.
In terms of the final residual, it seems that the larger ${\rm dim}\,{\cal Z}$ the better the performances. However, the gain might become limited as all the cases adopting the two-level preconditioner saturate on a similar level of residual $\|\hat{r}^{(k)}\| \sim 10^{-5} \div 10^{-6}$ . Moreover, the case with the largest basis, ${\rm dim}\,{\cal Z} = 50$, seems to enter a regime of high numerical noise that causes it to deliver worse final convergence compared to the other cases. Therefore, this numerical experiment suggests that, for fixed Arnoldi threshold, it is not recommended to strive for the largest deflation basis because an intermediate value can yield comparable (or even better) final convergence and consume less memory and CPU time – both proportional to ${\rm dim}\,{\cal Z}$.
In our set up, the case with $\epsilon_\lambda = 0.2$ (corresponding to ${\rm dim}\,{\cal Z} = 25$) delivers slightly more accurate estimate and will be our value of choice in the rest of this paper. As the threshold of $10^{-6}$ is commonly adopted in the CMB map-making procedures for the convergence, these residuals are already quite satisfactory. Moreover, they are also already nearly two orders of magnitude better than what can be achieved with the standard, block-diagonal preconditioner.

We would like to make sure that better precision could be reached if needed.  In the right panel of the figure we fix ${\rm dim}\,\mathcal{Z}$ to $25$ and show how the performances change as the  Arnoldi tolerance threshold decreases. The more we decrease the Arnoldi threshold, the lower value we get for the final residuals – for the reason discussed above, the first few tens of iterations are affected by the the fraction of eigenvectors retained rather than $\epsilon_{Arn}$ itself. Choosing $\epsilon_{Arn} \sim 10^{-6}$ seems already sufficient as it allows reaching residuals level as low as $10^{-8}$, but even lower residuals are reached by decreasing further  $\epsilon_{Arn}$. In particular, we did not reach any saturation when we let the Arnoldi converge completely, i.e. when $\epsilon_{Arn} = 10^{-14}$. We might be tempted to always use such a low threshold to build the preconditioner for our CG solver. However, when we push the Arnoldi to a given threshold we are basically solving the system to the same residual threshold with the GMRES algorithm. Therefore, building a two-level preconditioner for a given system using a value of $\epsilon_{Arn}$ much lower than the target CG residual is not meaningful.

We thus conclude that, in order to achieve a very accurate solution  (PCG residual tolerance $\sim 10^{-7}$ or better) by means of the two-level preconditioner, the Arnoldi algorithm has to converge within a tolerance of $\epsilon_{Arn}< 10^{-6}$, and ${\rm dim}\,{\cal Z}= 20\div 30$ eigenvectors are required to build the deflation subspace.

§.§ Exploitation of a precomputed two-level preconditioner

In the previous section, we have shown that building a two-level preconditioner with a fully converged Arnoldi algorithm gives the best CG convergence rate. Building such a preconditioner may not always be desirable for a single map-making run, given the extra numerical cost. Nonetheless, in this section we show that it typically not only leads to significant performance gains when many similar map-making runs are to be performed, but in process yields often better solutions for some single runs.

§.§.§ Divide-and-Conquer map-making of one season of observation

We now explore a different scenario, the so-called divide-and-conquer map-making, in which we solve for many map-making problems with a system matrix $A$ and right hand side (RHS)  $b$ that are similar but not equal.
CMB experiments can get the best possible map out of their observation only if they analyze the full data set at once. Nevertheless, splitting the full data volume into smaller groups and producing their maps independently can enormously reduce the computational complexity of map-making permitting to capitalize on the embarrassingly parallel character of this approach, while still producing high quality maps.

median        percentile

4*SP    $\mbd $     $ 3\cdot 10^{-5}$     0.0001    $8.1^{\dag}$

$\md, \epsilon_{Arn}=10^{-6}$     $9\cdot 10^{-6}$    $3 \cdot 10^{-5}$      $12.6^{\dag}$

$\md  $ simpl.      $7\cdot 10^{-7}$   $2\cdot 10^{-6}$      $10^{\dag}$

$\md,  \epsilon_{Arn}=10^{-12}$     $1\cdot 10^{-7}$   $3\cdot 10^{-6} $     $20.8^{\dag}$

2*LP    $\mbd $     $\si{\num{0.001}}$     0.0004	    $10.8$

$\md  $ simpl. 		    $3\cdot 10^{-7}$   $1\cdot 10^{-7}$     $8.5$

2*VLP    $\mbd $     $10^{-5}$     $ 10^{-5}	 $      $19.2$

$\md  $ simpl. 		     $ 10^{-7} $   $ 10^{-7}$      $15.9$

Median, $1\sigma$ statistics  of residual norms, $\parallel\hat{ r} ^{(k)} \parallel $ and computational cost of PCG runs for  different scanning strategies.We consider $p_{16th}$ and $p_{84th}$ respectively the  16-th and 84-th percentiles as $1\sigma$ upper and   lower bounds.  In the fourth column we quote $(p_{84th}-p_{16th})/2$. $^{\dag}$Values rescaled from Edison  to Cori  computational system  to a better compare performances.

Histograms of (left) residual norms and (right) iteration steps of PCG runs performed on simulated data on SP, LP and VLP respectively in top, center, bottom panel. Shown with  blue bars are the histogram related to PCG runs with $\mbd$, orange bars runs with $\md$  applied with the simplified approach.  In top panel are also shown PCG runs applied with the active approach with $\epsilon_{Arn}=10^{-6} \,(\,10^{-12}\, ) $ as red (green) bars.

In the context of the ground based experiments which typically scan the same sky area repetitively multiple times, these smaller map-making problems can be defined in such way that their system matrices $A$ have similar properties.

We explore the performances of the two-level preconditioner in this context starting from simulation of a two season data set of SP.
For this scanning strategy  each CES lasts about 15 minutes and  we split the whole observation into 250 subsets consisting of 27 CESs. This subgroup roughly corresponds to all the data taken in a given day.  Each processing element performs a PCG run on one of such subsets, which is characterized by $N_t\sim \si{\num{d8}}$ and $N_p \sim \si{\num{4d4}} $. Given these numbers, we can perform as many as two PCG runs per node of the Edison computing system[<http://www.nersc.gov/users/computational-systems/edison>], which provides 64 GB of memory.

We consider different types of two-level preconditioner runs

* The “Active” approach: the Ritz eigenpairs are computed for each subset of data. We use ${\rm dim}\,{\cal Z} = 25$ ($\epsilon_\lambda = 0.2$) but consider two values for $\epsilon_{Arn}$, $10^{-6}$ and $10^{-12}$. The former corresponds to the prescription we have given in Sect. <ref> for the single map-making run. The latter corresponds to the best preconditioner we can have with this technique. We explained earlier that it is not meaningful to build such a preconditioner for a PCG run, but it provides a useful limit case to compare against.

* The “Simplified” approach: the Ritz eigenpairs are computed only for one subgroup, using ${\rm dim}\,{\cal Z} = 25$ and $\epsilon_{Anr} = 10^{-12}$. The $\md $ built from this eigenvector basis is then applied to the rest of the whole data set. This approach is computationally cheaper than the active one, but it is less specific and it can work only in the case where the computed  deflation basis is very  representative of the whole dataset.

PCG error norms, $e^{\left(k\right)}$, as defined in Eq. (<ref>),   for a group of 27 signal-only CESs observing SP (left), for a group of 7 signal-only CESs observing LP (center), for a group of 3 signal-only CESs observing VLP (right).

Fig. <ref> reports the histograms of the PCG performances, the left and right panel respectively show the residual at the last iteration and the total number of iterations performed, the rightmost bin of the latter collects the runs that did not meet the convergence criterion of $10^{-6}$ within 100 iterations. As expected,  the degeneracies in the system matrix play a role in the runs involving $\mbd$ (blue bars), so that only $\sim 30\%$ converges below a threshold of $10^{-6}$ and most of the runs converges at higher residual tolerances. As far as $\md$ is concerned, the most remarkable result is that the simplified approach (orange histogram) perform considerably better than the active case with $\epsilon_{Arn} = 10^{-6}$ (red bars) and nearly as good as the one with  $\epsilon_{Arn} = 10^{-12}$ (green bars). This means that, rather than using a preconditioner tailored on the subset of the PCG run, it is important characterize well the most degenerate modes, even on a slightly different system: the deflation basis  computed from a subset of data is very well representative  for the  whole season data.
This demonstrates the simplified approach provides a very effective way of implementing the divide-and-conquer map-making run in the context of the ground-based observations.

Bin-by-bin comparison of power spectra differences defined in Eq. (<ref>) computed  from maps estimated at several iteration steps observing SP (top),  LP (center) and VLP (bottom). Dotted-blue and dashed-green lines refer respectively to  $E-$ and $B-$modes computed from maps reconstructed with Jacobi preconditioner.  Solid orange and red lines refer respectively to $E-$ and $B-$modes computed from maps reconstructed with $\md$. We choose different multipole bins to emphasize the convergence behavior of large, intermediate and small angular scales.

We further test this approach by applying it to observations covering a larger fraction of sky, the LP and VLP scanning strategies.  As summarized in table <ref> the noise level in the LP case  is about 4 times lower than the SP one: our aim is to probe the performances of this methodology in the perspective of the sensitivities that will be achieved by the forthcoming ground based experiments. For LP the length of one CES is larger than SP, each one lasting for 4 hours, and usually we simulate one or two CESs per day, depending on the seasonal availability of the patch above the horizon.  We obtain a data set consisting of 350 CESs and we chunk it into 50 subgroups made of 7 CESs.
Given the memory size of NERSC computing system  Cori ($128$ GB)[<http://www.nersc.gov/users/computational-systems/cori>]  we can run one chunk of data  made by 7 CES per node so that we distribute the seasonal data set across 50 nodes.
We construct the two-level preconditioner by taking one of these subsets and running the  Arnoldi iterations up to the numerical convergence,  $\epsilon_{Arn}=10^{-12}$. We retain the Ritz eigenvectors related to the eigenvalues smaller than $0.2$. This yields a deflation subspace with size $dim(\mathcal{Z}) =28$. We then apply the two-level preconditioner to the whole dataset.
The comparison of performances between the PCG run with $\mbd$ and $\md$ is shown in fig. <ref>, respectively with blue and orange bars.  Even in this case, we adopted the same tolerance,  $10^{-6}$,  and  we end the PCG iterations when this  tolerance  is not  achieved within 100 iterations. While most of the $\mbd$ runs do not converge, $\md$ runs converged within a median value of $ 44$ iterations.

PCG runs performed for 100 MC simulated data observing the VLP. Notice that all the 100 MC runs converged below a tolerance of $10^{-7}$ within $\sim  15$ iterations.

The last case to be analyzed is the VLP, which targets $20\%$ of the sky with a sensitivity typical of future CMB observatories.  Similarly to the SP and LP case, we compare the performances of the two preconditioners applied to one year of signal-plus-noise observations – a total of $300$ simulated 8-hour long CESs, grouped in 100 subsets of 3 CESs. We run the map-making solver on 100 processing elements distributed on 100 nodes of the Cori.
Consistently with the previous LP case, we apply the $\md$ with the simplified approach with an  Arnoldi tolerance of $\epsilon_{Arn}=10^{-12}$. In this case, the deflation subspace is spanned by 15 Ritz eigenvectors. As shown in Fig. <ref>, also in this case while the $\mbd$ rarely  converges to the fixed tolerance ($10^{-6}$),  whereas $\md$  converges within few tens of iterations.

Further details about the convergence statistics and total computational cost of SP, LP and VLP can be found in table <ref>.

§.§ Real space convergence

We analyze the convergence performances of the two-level and block-diagonal preconditioners using the norm of the difference as defined in Eq. (<ref>). Compared to the standard PCG residuals, this metric emphasizes more the eigenvectors of the system matrix with low eigenvalues. As mentioned earlier this analysis requires the knowledge of the exact solution of the system. For this reason, we perform  signal-only simulations for a  subgroup of all  the  observational patches discussed in sect. <ref> and compare the performance of $\mbd$ and $\md$, computed with the  active approach.

The results, shown in Fig. <ref>, show that the two-level preconditioner is able to recover the solution to  some order of magnitude better precision with respect to  the one computed with the block-diagonal methodology. The fact that the latter saturates very quickly at a value way higher then the PCG residual emphasizes further that its convergence is hindered by the nearly degenerate modes, which are downweighted in the PGC residuals shown in the other plots,  Fig <ref>. Moreover, the fact that the  saturation levels differs case by case  in Fig. <ref>, could be  due to the presence of different degeneracies depending on the considered observational patch.

§.§ Convergence at the power spectrum level

We investigate a scale-dependence of the reconstructions by analyzing the signal-only study cases considered in the previous section and perform the bin-by-bin power spectra comparison of the residuals as shown in Fig <ref>.

For the  SP case, the two-level preconditioner converges to the threshold of $10^{-7}$ within 40 iterations, whereas the case with $\mbd$  do not within 100 iterations, i.e., the maximum allowed in these runs. We consider the bins that are usually considered in the analysis of patches as small as the SP.  As one can notice from Fig. <ref> (top), the  solution computed with $\mbd $ encodes an extra-bias which is order of few percentages sub-dominant with respect to the variance of the signal itself, meaning that the quality of the map reconstructed with the $\mbd $ is acceptable as far as small angular scales are concerned. Moreover, this is somewhat expected since the larger angular scales are not constrained by the $\mbd$ and are the responsible of the long mode plateau we described in sect. <ref>.  Those scales are anyway unconstrained due to the small sizes of the patch.

LP allows us to probe larger scales, where the primordial gravitational wave B-mode signal is expected to peak. The solution computed with $\mbd$ (which does not converge within 100 iterations) shows a $\sim 10 \%$ bias at the largest angular scales ( in the first two bins, namely $\ell_b = 50-100, 150-200$ in Fig. <ref>(center panel)), whereas the bias is  not present into  the solution computed with $\md$. This result becomes even more remarkable given that at these scales the signal is likely to be dominated by foreground emission, therefore the same fractional bias in the power-spectrum can be comparable with the whole signal from primordial $B$-modes.  In terms of the norm of the standard residuals these results demonstrate that high precision convergence needs to be attained in order to ensure a sufficient precision of the recovered sky signal on all, and specifically on the largest accessible, angular scales.

We observe a similar behavior with the spectra computed for VLP, Fig.<ref>(bottom). In particular, we focus onto  large scales since the size of the patch is big enough to probe the angular scales related to the reionization peak of both $E$- and $B$-modes. We notice that the first two multipole bins $\ell_{b}=0-50$ and $\ell_{b}=50-100$ are reconstructed up to percentage level with the two-level preconditioner, whereas the power spectra computed with the block-diagonal one contains a bias which may fluctuate between tens and few percentages.  The degree and subdegree angular scales are similarly reconstructed as in the LP case.

§.§.§ Monte Carlo simulations

All modern CMB experiments produce or validate their statistical and systematic uncertainties using a large number of simulations. Typically, each of them solves for a map-making system that has the same system matrix $A$ but different RHS $b$ (i.e., the same scanning strategy and data processing but different synthetic time stream).
We consider an observation composed of 3 CESs covering the VLP. We produce 100 Monte Carlo (MC) with not only uncorrelated noise, but even a signal generated using different random seeds from the same CAMB power spectra. We take one of these simulations and build a two-level preconditioner from a fully converged Arnoldi run. We then apply the same preconditioner to all the simulations. As shown in Fig. <ref>, all the 100 MC runs converged to a residual tolerance $< 10^{-7}$ within $\sim 15$ iterations and with a staggering narrow dispersion. This result on one hand shows how powerful a two-level preconditioner can be when MC simulations are to be performed, on the other it means that the degeneracies preventing the convergence with the standard preconditioner are not due to the signal or the presence of noise, but mostly due to the scanning strategy and the filtering applied to the time stream.

§ SUMMARY AND CONCLUSIONS

In this work we described an implementation of a novel class of iterative solvers, the two-level preconditioners,  $\md$, [Grigori et al., 2012, Szydlarski et al., 2014] in the context of the map-making procedure applied to data sets filtered at the time domain level. We have discussed the details of the construction of the new preconditioner and proposed a simplified, ”divide and conquer”, embarrassingly parallel implementation of the method, which can be adequate for an analysis of current and future, ground-based observations.
We have tested this new implementation of this novel methodology on three different simulated data sets in the cases when filtering operators typical of the ground experiments, have been applied. We have compared the performance of the method with that of the standard PCG solver based on the Jacobi preconditioner.

We have found that in all the studied cases, the two-level preconditioner, $\md$, have performed better both in terms of the attained precision and the number of required iterations, allowing typically reaching the residuals on order $10^{-7}$ within $20\div 40$ iterations. The standard approach yields residuals an order of magnitude or more higher within as many as $100$ iterations. We show that reaching such high precision of the reconstructed maps is required in order to constrain the  large angular scales of the B-mode polarization. Indeed, the new approach consistently produces maps typically within $20\div40$ iterations, which display negligible reconstruction bias of all and in particular the longest modes as represented in the maps. In the contrary, the maps derived with the standard solver with the maximal number of iterations set to $100$ show typically a $1-20 \%$ bias at all scales.

We thus conclude that producing highly accurate maps of the polarized CMB anisotropies from the filtered data of the ground-based experiments may call for more advanced iterative solvers than the standard PCG solver with the Jacobi preconditioner. The presented here, two-level preconditoner offers significantly better performance and could be a method of choice for such applications in the future. These advantages come however at the additional cost needed to construct the preconditioner. Therefore this method can be particularly useful in the cases of large MC simulations, where the additional cost is offset by the solver's superior performance.

We acknowledge use of camb, healpix, s4cmb, and x2pure software packages. We thank Julien Peloton for his help with simulations, and Josquin Errard and Maude Le Jeune for helpful comments and suggestions.
RS acknowledges support of the French National Research Agency (ANR) contract ANR-17-C23-0002-01 (project B3DCMB).

This research used resources of the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. GF acknowledges the support of the CNES postdoctoral program. CB, DP acknowledge support from the RADIOFOREGROUNDS project, funded by the European Commission’s H2020 Research Infrastructures under the Grant Agreement 687312, and the COSMOS Network from the
Italian Space Agency; CB acknowledges support by the INDARK INFN Initiative.

[Abazajian et al., 2016]
Abazajian, K. N., Adshead, P., Ahmed, Z., et al. 2016, ArXiv e-prints

[Ade et al., 2015]
Ade, P. A. R., Aghanim, N., Ahmed, Z., et al. 2015, Phys. Rev. Lett., 114,

[Arnold et al., 2014]
Arnold, K., Stebor, N., Ade, P. A. R., et al. 2014, in , Vol.
9153, Millimeter, Submillimeter, and Far-Infrared Detectors and
Instrumentation for Astronomy VII, 91531F

[Ashdown et al., 2009]
Ashdown, M. A. J., Baccigalupi, C., Bartlett, J. G., et al. 2009, ,
493, 753

[Benson et al., 2014]
Benson, B. A., Ade, P. A. R., Ahmed, Z., et al. 2014, in ,
Vol. 9153, Millimeter, Submillimeter, and Far-Infrared Detectors and
Instrumentation for Astronomy VII, 91531P

[BICEP2 Collaboration et al., 2016]
BICEP2 Collaboration, Keck Array Collaboration, Ade, P. A. R., et al.
2016, , 833, 228

[Borrill, 1999]
Borrill, J. 1999, ArXiv Astrophysics e-prints [astro-ph/9911389]

[Bunn et al., 2003]
Bunn, E. F., Zaldarriaga, M., Tegmark, M., & de Oliveira-Costa, A.
2003, , 67, 023501

[Cantalupo et al., 2010]
Cantalupo, C. M., Borrill, J. D., Jaffe, A. H., Kisner, T. S., &
Stompor, R. 2010, , 187, 212

[de Gasperis et al., 2005]
de Gasperis, G., Balbi, A., Cabella, P., Natoli, P., & Vittorio, N.
2005, , 436, 1159

[Doré et al., 2001]
Doré, O., Teyssier, R., Bouchet, F. R., Vibert, D., & Prunet, S.
2001, , 374, 358

[Golub & Van Loan, 1996]
Golub, G. H. & Van Loan, C. F. 1996, Matrix Computations (3rd Ed.) (Baltimore,
MD, USA: Johns Hopkins University Press)

[Górski et al., 2005]
Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, , 622,

[Grain et al., 2009]
Grain, J., Tristram, M., & Stompor, R. 2009, , 79, 123515

[Grigori et al., 2012]
Grigori, L., Stompor, R., & Szydlarski, M. 2012, in Proceedings of the
International Conference on High Performance Computing, Networking, Storage
and Analysis, SC '12 (Los Alamitos, CA, USA: IEEE Computer Society Press),

[Hanson et al., 2013]
Hanson, D., Hoover, S., Crites, A., et al. 2013, Physical Review
Letters, 111, 141301

[Henderson et al., 2016]
Henderson, S. W., Allison, R., Austermann, J., et al. 2016, Journal of
Low Temperature Physics, 184, 772

[Hivon et al., 2002]
Hivon, E., Górski, K. M., Netterfield, C. B., et al. 2002, ,
567, 2

[Huffenberger & Næss, 2018]
Huffenberger, K. M. & Næss, S. K. 2018, , 852, 92

[Keisler et al., 2015]
Keisler, R., Hoover, S., Harrington, N., et al. 2015, , 807, 151

[Lewis et al., 2000]
Lewis, A., Challinor, A., & Lasenby, A. 2000, Astrophys. J., 538, 473

[Lewis et al., 2001]
Lewis, A., Challinor, A., & Turok, N. 2001, Phys. Rev. D, 65, 023505

[Louis et al., 2017]
Louis, T., Grace, E., Hasselfield, M., et al. 2017, , 6, 031

[Naess et al., 2014]
Naess, S., Hasselfield, M., McMahon, J., et al. 2014, , 10, 007

[Oh et al., 1999]
Oh, S. P., Spergel, D. N., & Hinshaw, G. 1999, , 510, 551

[Planck Collaboration et al., 2016]
Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, ,
594, A13

[Poletti et al., 2017]
Poletti, D., Fabbian, G., Le Jeune, M., et al. 2017, , 600, A60

[Smith & Zaldarriaga, 2007]
Smith, K. M. & Zaldarriaga, M. 2007, Phys. Rev. D, 76, 043001

[Stompor et al., 2002]
Stompor, R., Balbi, A., Borrill, J. D., et al. 2002, , 65, 022003

[Szydlarski et al., 2014]
Szydlarski, M., Grigori, L., & Stompor, R. 2014, , 572, A39

[Tegmark, 1997]
Tegmark, M. 1997a, , 56, 4514

[Tegmark, 1997]
Tegmark, M. 1997b, , 480, L87

[The Polarbear Collaboration et al., 2017]
The Polarbear Collaboration, Ade, P. A. R., Aguilar, M., et al. 2017,
, 848, 121

[The Polarbear Collaboration et al., 2014]
The Polarbear Collaboration, Ade, P. A. R., Akiba, Y., et al.
2014a, Physical Review Letters, 112, 131302

[The Polarbear Collaboration et al., 2014]
The Polarbear Collaboration, Ade, P. A. R., Akiba, Y., et al.
2014b, , 794, 171

[Wright, 1996]
Wright, E. L. 1996, ArXiv Astrophysics e-prints [astro-ph/9612006]

§ THE ARNOLDI ALGORITHM

The Krylov subspace algorithms are based on the construction of a sequence of vectors naturally produced by the power method,  a class of  those, namely the  Minimal Residuals (MINRES) and the Generalized Minimal Residual (GMRES) methods [Golub & Van Loan, 1996], rely on the Arnoldi algorithm.
Our goal is to find an approximation to the eigenvalues of a matrix $B$ of  a generic linear system with RHS $b$:

\begin{equation}
\label{eq:newlinsys}
\end{equation}

The Arnoldi algorithm is an algorithm aimed at solving linear systems by projecting the system matrix onto a convenient Krylov subspace generated by the first $m$ vectors

\begin{equation}
\mathcal{K}_m (B, b)= span \lbrace b, Bb,B^2 b, \dots, B^{m-1} b \rbrace.
\label{eq:krylov}
\end{equation}

The  major steps of the algorithm are summarized in Algorithm <ref>.
Basic steps of the Arnoldi algorithm
[1]
: $r_0,\, w_1= r_0/\parallel r_0 \parallel$
j = 1$\rightarrow $ m
$v_j=B w_j$
i = 1$\rightarrow $ j
$h_{i,j}= (v_j , w_i)$
$v_j = v_j - \sum_{i=1} ^j h_{i,j} w_i$
$h_{j+1,j} = \parallel  v_j\parallel$
$h_{j+1,j}\leq \mathrm{tol}_{Arn}$
$w_{j+1} = v_j / h_{j+1,j}$

Hence, the output of the Arnoldi algorithm is an orthonormal basis $W^{(m)}= (w_1| w_2|\dots| w_m)$ (called the Arnoldi vectors ), together with a set of scalars $h_{i,j}$ (with $i,j=1,\dots,m$ and $i\leq j+1$) plus an extra-coefficient $h_{m+1,m}$. The former set of coefficients are the elements of an upper Hessenberg matrix $H_m $ with non-negative subdiagonal elements and is commonly referred as a m-step Arnoldi Factorization of $B$. If $B$ is Hermitian then $H_m$ is symmetric, real  and  tridiagonal and the vectors (columns of $W^{(m)}$) of the Arnoldi basis  are called Lanczos vectors. $B$ and $H_m$ are intimately related via:

\begin{equation}
BW^{(m)} = W^{(m)} H _m + h_{m+1,m} w_{m+1} e^\dag  _m,
\label{eq:factoriz}
\end{equation}

where $e_m$ is a $1\times m $ unit vector with 1 on the $m$-th component.

In other words, $H_m$ is the projection of $B$ onto the subspace generated by the Arnoldi basis $W^{(m)}$  within an error  given by $\tilde{W}_m= h_{m+1,m} w_{m+1} e^\dag  _m$.  The  iteration loop ends when this error term gets smaller than a  certain threshold $\epsilon_{Arn}$.

Using Eq. (<ref>), we can connect the eigenpairs of $B$ to the ones of $H_m$.
Let us consider an eigenpair of $H_m, \,(\lambda_i,y_i)$

H_m y_i = λ_i y_i.

The vector $v_i= W^{(m)} y_i$ then satisfies

\begin{equation}
\parallel B v_i - \lambda_i v_i \parallel= \parallel ( B W^{(m)} - W^{(m)} H_m ) v_i \parallel=  \parallel \tilde{W}_m v_i \parallel.
\label{eq:ritz}
\end{equation}

The eigenpairs of $H_m$ are therefore approximations of the eigenpairs of $B$ within an  error given by $\tilde{W}_{m+1}$. They are the so called Ritz eigenpairs  and they are very easy to compute since the size of $H_m$ is  $\lesssim \mathcal{O}(100)$. For the dataset considered in this work this is indeed the order of Arnoldi iterations required to reach a tolerance $\epsilon_{Arn}\sim 10^{-6}$ as it can be seen in Fig. <ref>. A typical distribution of the amplitude of the Ritz eigenvalues for different values of $\epsilon_{Arn}$ is shown in in Fig. <ref>.

The convergence residuals of the Arnoldi algorithm. Notice that after about $175 $ iterations the algorithm numerically converges.

PCG residuals for SP intensity maps different choices of $\epsilon_{Arn}$  and at a given $\mathrm{dim} \mathcal{Z}= 25$.

§ SOLVING FOR INTENSITY MAPS

Similarly to what we have done in Sect. <ref>, we further tested the two-level preconditioner on SP intensity-only maps. As it is shown in Fig. <ref>, in this case the time stream has to be filtered with  a higher  polynomial basis due to a larger  $f_{knee}$. We therefore filter the time streams up to the third order Legendre polynomials.

Histograms of Ritz eigenvalues of the matrix $\mbd A$ computed for several choices of Arnoldi tolerance, $\epsilon_{Arn}$.

Fig. <ref> shows the PCG residuals for different choices of $\epsilon_{Arn}$ and 	one can easily notice the similarity to the right panel of Fig. <ref>. This further indicates that our results are stable even when a more aggressive filter is applied to the data. Moreover, by looking at the blue-dashed line in  Fig. <ref>, the $\mbd$  residuals saturate at a higher threshold with respect to the polarization case (Fig.<ref>), remarking  the presence  of different degeneracies present when intensity maps are involved. However, the two-level preconditioner does not suffer of this effect and once the Ritz eigenvector basis is very well approximated, by running the Arnoldi algorithm to tolerances below $10^{-6}$,   it converges to $10^{-7}$ within $\sim 40$ iterations.


# indicnlp@kgp at DravidianLangTech-EACL2021: Offensive Language
Identification in Dravidian Languages

Kushal Kedia
IIT Kharagpur
<EMAIL_ADDRESS>
Abhilash Nandy
IIT Kharagpur
<EMAIL_ADDRESS>

###### Abstract

The paper presents the submission of the team indicnlp@kgp to the EACL 2021
shared task “Offensive Language Identification in Dravidian Languages”. The
task aimed to classify different offensive content types in 3 code-mixed
Dravidian language datasets. The work leverages existing state of the art
approaches in text classification by incorporating additional data and
transfer learning on pre-trained models. Our final submission is an ensemble
of an AWD-LSTM based model along with 2 different transformer model
architectures based on BERT and RoBERTa. We achieved weighted-average F1
scores of $0.97$, $0.77$, and $0.72$ in the Malayalam-English, Tamil-English,
and Kannada-English datasets ranking ${1}^{st}$, ${2}^{nd}$, and $3^{rd}$ on
the respective tasks.

##  1 Introduction

Offensive language identification is a natural language processing (NLP) text
classification task where the goal is to moderate and reduce objectionable
social media content. There has been a rapid growth in offensive content on
social media and the number of users from different ethnicities and cultures
worldwide. A significant portion of offensive content is specifically targeted
at various individuals and minority & ethnic groups. Consequently, the
identification and classification of these different kinds of foul language
are receiving increased importance. Dravidian languages like Kannada,
Malayalam, and Tamil Raja and Chakravarthi (2020) are low-resourced making
this task challenging. Training embeddings of words has previously been a
common approach employed in text classification tasks. However, transfer
learning approaches in deep learning Mou et al. (2016) have been shown
unsuccessful or requiring extremely large collections of in-domain documents
to produce strong results Dai and Le (2015).

Further, in a multilingual culture, code-mixing is a prevalent practice and
code-mixed texts are sometimes written in native scripts. Due to the
complexity introduced by code-switching at multiple linguistic levels, systems
trained on monolingual data can fail on code-mixed data. While multilingual
versions of transformer models have been shown to perform remarkably well,
even in zero-shot settings Pires et al. (2019), a zero-shot transfer may
perform poorly or fail altogether Søgaard et al. (2018). This is when the
target language, here code-mixed Dravidian data, is different from the source
language, mainly monolingual. In our work, we tackle these problems by
exploiting additional datasets for fine-tuning our models and using effective
transfer learning techniques. Our code and experiments are available on
GitHub111https://github.com/kushal2000/Dravidian-Offensive-Language-
Identification for generating reproducible results on our models.

##  2 Task Description and Datasets

The purpose of this task is to classify offensive language material gathered
from social media from the a set of code-mixed posts in Dravidian Languages.
The systems have to classify each post into one of the 6 labels:

  * •

not offensive

  * •

untargeted offense

  * •

offense targeted at an individual

  * •

offense targeted at a group

  * •

offense targeted at someone else

  * •

not in intended language

There is also a significant class imbalance in all the datasets representing a
real-world situation. This shared task presents a new gold standard corpus for
offensive language identification of code-mixed text in three Dravidian
languages: Tamil-English Chakravarthi et al. (2020b), Malayalam-English
Chakravarthi et al. (2020a), and Kannada-English Hande et al. (2020). The
Malayam dataset does not contain the offense targeted at someone else tag. The
posts can contain more than one sentence, but the average number of sentences
is $1$. The Tamil and Malayalam datasets are considerably large containing
over $30k$ and $20k$ annotated comments while the Kannada dataset is
relatively smaller with almost $8k$ annotations. Apart from the dataset
supplied by the organizers, we also use a monolingual English Offensive
Language Identification Dataset (OLID) Zampieri et al. (2019a) used in the
SemEval-2019 Task 6 (OffensEval) Zampieri et al. (2019b). The dataset contains
the same labels as our task datasets with the exception of the not in intended
language label. The one-to-one mapping between the labels in OLID and it’s
large size of $14k$ tweets makes it suitable for aiding the transfer learning
detailed in Section 3.3.

##  3 Methods

A variety of methods are experimented on the datasets to provide a complete
baseline. In Section 3.1, we describe our implementation of three traditional
machine learning classifiers; Multinomial Naive Bayes, Linear Support Vector
Machines (SVM) and Random Forests. These approaches work well on small
datasets and are more computationally efficient than deep neural networks.
Their performance is similar to the models described in the latter sections in
the absence of pretraining and additional data. In Section 3.2, our Recurrent
Neural Network (RNN) models are explained. We have compared an LSTM model
using word-level embeddings trained from scratch, to an ULMFiT model, an
effective transfer learning approach for language models. Finally, in Section
3.3, we discuss transformer architectures using their cross-lingual pretrained
models which can be data intensive during fine-tuning but provide the
strongest results on our datasets.

###  3.1 Machine Learning Classifiers

Dataset Preprocessing The datasets are preprocessed by removing punctuation,
removing English stop words, removing emojis, and lemmatizing the English
Words. The Natural Language Toolkit library Bird and Loper (2004) was used for
lemmatization and removing stop words. A word vocabulary is made and
vocabulary-length vectors containing counts of each word are used to represent
each individual input. Based on the Mutual Information scores of each word,
feature selection is done to reduce the vocabulary size.

Hyperparameters For all three models, the number of words selected using the
top Mutual Information scores was varied from $1000$ to the length of the
vocabulary. Further hyperparameters were specific to the SVM and Random
Forest. The random state, the regularisation parameters, and max iterations
were tuned for the SVM and the number of decision trees used was the only
hyperparameter in the case of random forests.

###  3.2 RNN Models

Vanilla LSTM To set a baseline for an RNN approach, we build word embeddings
from scratch using just the individual datasets. For this, we selected the top
$32,000$ occurring words in each dataset for one-hot encoding, which is passed
through an embedding layer to form $100$-dimension word vectors. A spatial
dropout of $0.2$ followed by a single LSTM cell and a final softmax activation
forms the rest of the model. While the results for larger datasets are
marginally better than the previous section, they are worse in comparison to
the transfer learning approach.

ULMFiT Transfer learning has been shown to perform well in text classification
tasks. Usually, language models are trained on large corpora and their first
layer, i.e, the word embeddings are fine-tuned on specific tasks. This
approach has been a very successful deep learning approach in many state of
the art models. Mikolov et al. (2013) However, Howard and Ruder, 2018 argue
that we should be able to do better than randomly initializing the remaining
parameters of our models and propose ULMFiT: Universal Language Model Fine-
tuning for Text Classification. For the Dravidian languages in this task, the
problem of in-domain data collection for effective transfer is also
significant especially in the domains of hate speech. ULMFiT provides a robust
framework for building language models from moderate corpora and fine-tunes
them on our specific tasks.

Language Models & Corpora We make use of language models open-soruced by the
team gauravarora Arora (2020) in the shared task at HASOC-Dravidian-CodeMix
FIRE-2020 Mandl et al. (2020). They build their corpora for language modelling
from large sets of Wikipedia articles. For Tamil & Malayalam languages, they
also generate code-mixed corpora by obtaining parallel sets of native,
transliterated and translated articles and sampling sentences using a Markov
process, which has transition probabilities to $3$ states; native, translated
and transliterated. For Kannada, only a native script corpus is available and
we had to transliterate our code-mixed dataset to Kannada to match their
language model. The models are based on the Fastai Howard and Gugger (2020)
implementation of ULMFiT. Pre-trained tokenizers and language models are
available on Github. 222github.com/goru001/nlp-for-tanglish
333github.com/goru001/nlp-for-manglish 444github.com/goru001/nlp-for-kannada

Preprocessing & Model Details Basic preprocessing steps included lower-casing,
removing punctuations and mentions, etc. Subword tokenization using unigram
segmentation is implemented which is reasonably resilient to variations in
script and spelling. The tokenization model used is
SentencePiece555github.com/google/sentencepiece. The language model is based
on an AWD-LSTM Merity et al. (2018), a regular LSTM cell with additional
parameters related to dropout within the cell. The text classification model
additionally uses two linear layers followed by a softmax on top of the
language model. To tackle the difference in distributions of the target
datasets and the pretraining corpora, ULMFiT proposes using 1) discriminative
fine-tuning, i.e, layers closer to the last layer have higher learning rates,
2) slanted triangular learning rates which increase aggressively during the
start of training and then decay gradually and 3) gradual unfreezing, i.e,
instead of learning all layers of the model at once, they are gradually
unfrozen starting from the last layer. The combination of these techniques
leads to robust transfer learning on our datasets.

###  3.3 Transformer Models

In recent years, transformer networks like the Bidirectional Encoder
Representation from Transformer (BERT) Devlin et al. (2019) and its variant
RoBERTa Liu et al. (2019) have been used successfully in many offensive
language identification tasks. For our work, we use the already pre-trained
cross-lingual versions of these models available in the
HuggingFace666https://huggingface.co/ library. Specifically, we use the bert-
base-multilingual-cased model, mBERT trained on cased text in $104$ languages
from large Wikipedia articles, and the xlm-roberta-base model, XLM-R Conneau
et al. (2020) trained on $100$ languages, using more than two terabytes of
filtered CommonCrawl data. Both of these models were originally trained on a
masked language modelling objective and we fine-tune them on our specific
downstream text classification tasks.

Transfer Learning The core principle of the transfer learning approach is to
use a pretrained transformer model for training a classification model on a
resource-rich language first, usually English, and transfer the model
parameters on a less resource-rich language. For this approach, we concatenate
all 3 code-mixed datasets as well as the OLID dataset. Our results do not
change significantly on transliteration of all datasets to Roman script. The
not in intended language label is also removed for fine-tuning on the combined
dataset since this label does not represent the same meaning across the
datasets. We then use these learned model weights replacing the final linear
layer to include the additonal removed label not in intended language. This
kind of fine-tuning approach has been shown to increase the performance of
various scarce-resourced languages such as Hindi and Bengali, etc. Ranasinghe
and Zampieri (2020).

Figure 1: Transformer Model Architecture

Model Architecture We restrict the maximum length of the input sentences to be
$256$ by truncation and zero-padding. As shown in Fig 1, using the contextual
embeddings from the last hidden states of all tokens in the sentence, we build
a vector representation by concatenating the max pooling and mean pooling of
these hidden states. Correspondingly, the dimension of the final sentence
representation is $1536$ x $1$ . This is passed through a linear layer with a
dropout of $0.3$. The learning rate for all fine-tuning was fixed as $2e^{-5}$
and batch size was $32$. The only preprocessing step before feeding the input
to the transformer tokenizers was the replacement of emojis by their
description in English. This is done because the tokenizers might not
recognize the emojis, but they contain useful information about the sentiment
of the sentence.

##  4 Experiments and Results

The assessment metric for this task is Weighted-F1, which is the F1 score
weighted by number of samples in all the classes. The datasets are divided
into train, validation, and test sets in an approximately 8:1:1 ratio. The
test labels are hidden and only available to us after the evaluation phase is
over. We strictly train our models on the train set using the scores on the
validation set for hyperparameter tuning. We have reported the results of our
various models on the validation set. The average-ensemble of our top three
performing models is submitted finally and we also report it’s scores on the
validation and test set - using the scores we are ranked with in the task
leader board.

Model | T | M | K
---|---|---|---
Random Forest | 0.69 | 0.94 | 0.62
Naive Bayes | 0.74 | 0.94 | 0.64
Linear SVM | 0.74 | 0.95 | 0.65
Table 1:  Weighted-F1 scores for ML models on Tamil (T), Malayalam (M) and
Kannada (K) datasets.

Table 1 showcases the scores we have obtained for standard machine learning
algorithms. Out of the three traditional machine learning algorithms, the
Linear SVM model is best across all three datasets. The results in Table 2
summarize our RNN approaches where ULMFiT is markedly superior. The
performance of our transformer models detailed in Table 3 considers two
settings, one without transfer learning and one with transfer learning using
the OLID and other Dravidian code-mixed datasets in conjunction.

Model | T | M | K
---|---|---|---
Vanilla LSTM | 0.74 | 0.95 | 0.64
ULMFiT | 0.76 | 0.96 | 0.71
Table 2:  Weighted-F1 scores for RNN models on Tamil (T), Malayalam (M) and Kannada (K) datasets. Model | T | M | K
---|---|---|---
mBERT | 0.74 | 0.95 | 0.66
XLM-R | 0.76 | 0.96 | 0.67
mBERT (TL) | 0.75 | 0.97 | 0.71
XLM-R (TL) | 0.78 | 0.97 | 0.72
Table 3:  Weighted-F1 scores for transformers on Tamil (T), Malayalam (M) and
Kannada (K) datasets. TL indicates transfer learning using OLID and other
datasets.

The results on the validation set of our transfer-learnt XLM-R model are the
best across all $3$ datasets and is followed closely by the transfer learnt
multilingual BERT model and the ULMFiT model. We finally submit an average
ensemble of these three models and our results on the validation set and the
test set used in the final task evaluation are also enlisted in Table 4 below.

Model | T | M | K
---|---|---|---
avg-Ensemble (V) | 0.78 | 0.97 | 0.73
avg-Ensemble (T) | 0.77 | 0.97 | 0.72
Table 4:  Final weighted-F1 scores using average ensembling on Tamil (T),
Malayalam (M) and Kannada (K) validation (V) and test (T) datasets.

##  5 Conclusion

This paper describes various approaches for offensive language identification
in three code-mixed English-Dravidian language datasets. We also discuss the
final system submitted by the indicnlp@kgp team, which ranks first, second,
and third on the competition’s three tasks. The benefit of pre-trained
language models was shown by the significant improvement in results using a
robust transfer learning framework (ULMFiT) compared to a vanilla LSTM model
trained from scratch. Transformer networks’ performance also improved when all
the Dravidian language datasets were combined. This suggests that learning
from one Dravidian language may help in zero-shot or few-shot transfer to
other new Dravidian languages. In future works, we wish to explore these
effects in more detail.

## References

  * Arora (2020) Gaurav Arora. 2020.  Gauravarora@hasoc-dravidian-codemix-fire2020:pre-training ulmfit on synthetically generated code-mixed data for hate speech detection.  In _FIRE-2020 (Working Notes)_. CEUR, Hyderabad, India.
  * Bird and Loper (2004) Steven Bird and Edward Loper. 2004.  NLTK: The natural language toolkit.  In _Proceedings of the ACL Interactive Poster and Demonstration Sessions_ , pages 214–217, Barcelona, Spain. Association for Computational Linguistics.
  * Chakravarthi et al. (2020a) Bharathi Raja Chakravarthi, Navya Jose, Shardul Suryawanshi, Elizabeth Sherly, and John Philip McCrae. 2020a.  A sentiment analysis dataset for code-mixed Malayalam-English.  In _Proceedings of the 1st Joint Workshop on Spoken Language Technologies for Under-resourced languages (SLTU) and Collaboration and Computing for Under-Resourced Languages (CCURL)_ , pages 177–184, Marseille, France. European Language Resources association.
  * Chakravarthi et al. (2020b) Bharathi Raja Chakravarthi, Vigneshwaran Muralidaran, Ruba Priyadharshini, and John Philip McCrae. 2020b.  Corpus creation for sentiment analysis in code-mixed Tamil-English text.  In _Proceedings of the 1st Joint Workshop on Spoken Language Technologies for Under-resourced languages (SLTU) and Collaboration and Computing for Under-Resourced Languages (CCURL)_ , pages 202–210, Marseille, France. European Language Resources association.
  * Conneau et al. (2020) Alexis Conneau, Kartikay Khandelwal, Naman Goyal, Vishrav Chaudhary, Guillaume Wenzek, Francisco Guzmán, Edouard Grave, Myle Ott, Luke Zettlemoyer, and Veselin Stoyanov. 2020.  _Unsupervised Cross-lingual Representation Learning at Scale_.  Association for Computational Linguistics.
  * Dai and Le (2015) Andrew M Dai and Quoc V Le. 2015.  Semi-supervised sequence learning.  In _Advances in Neural Information Processing Systems_ , volume 28, pages 3079–3087. Curran Associates, Inc.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 4171–4186, Minneapolis, Minnesota. Association for Computational Linguistics.
  * Hande et al. (2020) Adeep Hande, Ruba Priyadharshini, and Bharathi Raja Chakravarthi. 2020.  KanCMD: Kannada CodeMixed dataset for sentiment analysis and offensive language detection.  In _Proceedings of the Third Workshop on Computational Modeling of People’s Opinions, Personality, and Emotion’s in Social Media_ , pages 54–63, Barcelona, Spain (Online). Association for Computational Linguistics.
  * Howard and Gugger (2020) Jeremy Howard and Sylvain Gugger. 2020.  Fastai: A layered api for deep learning.  _Information_ , 11(2):108.
  * Howard and Ruder (2018) Jeremy Howard and Sebastian Ruder. 2018.  Universal language model fine-tuning for text classification.  In _Proceedings of the 56th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 328–339, Melbourne, Australia. Association for Computational Linguistics.
  * Liu et al. (2019) Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  Roberta: A robustly optimized bert pretraining approach.
  * Mandl et al. (2020) Thomas Mandl, Sandip Modha, Anand Kumar M, and Bharathi Raja Chakravarthi. 2020\.  Overview of the hasoc track at fire 2020: Hate speech and offensive language identification in tamil, malayalam, hindi, english and german.  In _Forum for Information Retrieval Evaluation_ , FIRE 2020, page 29–32, New York, NY, USA. Association for Computing Machinery.
  * Merity et al. (2018) Stephen Merity, Nitish Shirish Keskar, and Richard Socher. 2018.  Regularizing and optimizing LSTM language models.  In _International Conference on Learning Representations_.
  * Mikolov et al. (2013) Tomas Mikolov, Ilya Sutskever, Kai Chen, Greg S Corrado, and Jeff Dean. 2013.  Distributed representations of words and phrases and their compositionality.  In _Advances in Neural Information Processing Systems_ , volume 26, pages 3111–3119. Curran Associates, Inc.
  * Mou et al. (2016) Lili Mou, Zhao Meng, Rui Yan, Ge Li, Yan Xu, Lu Zhang, and Zhi Jin. 2016.  How transferable are neural networks in NLP applications? In _Proceedings of the 2016 Conference on Empirical Methods in Natural Language Processing_ , pages 479–489, Austin, Texas. Association for Computational Linguistics.
  * Pires et al. (2019) Telmo Pires, Eva Schlinger, and Dan Garrette. 2019.  How multilingual is multilingual BERT? In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , pages 4996–5001, Florence, Italy. Association for Computational Linguistics.
  * Raja and Chakravarthi (2020) Bharathi Raja and Asoka Chakravarthi. 2020.  _Leveraging orthographic information to improve machine translation of under-resourced languages_.  NUI Galway.
  * Ranasinghe and Zampieri (2020) Tharindu Ranasinghe and Marcos Zampieri. 2020.  Multilingual offensive language identification with cross-lingual embeddings.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 5838–5844, Online. Association for Computational Linguistics.
  * Søgaard et al. (2018) Anders Søgaard, Sebastian Ruder, and Ivan Vulić. 2018.  On the limitations of unsupervised bilingual dictionary induction.  In _Proceedings of the 56th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 778–788, Melbourne, Australia. Association for Computational Linguistics.
  * Zampieri et al. (2019a) Marcos Zampieri, Shervin Malmasi, Preslav Nakov, Sara Rosenthal, Noura Farra, and Ritesh Kumar. 2019a.  Predicting the type and target of offensive posts in social media.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 1415–1420, Minneapolis, Minnesota. Association for Computational Linguistics.
  * Zampieri et al. (2019b) Marcos Zampieri, Shervin Malmasi, Preslav Nakov, Sara Rosenthal, Noura Farra, and Ritesh Kumar. 2019b.  SemEval-2019 task 6: Identifying and categorizing offensive language in social media (OffensEval).  In _Proceedings of the 13th International Workshop on Semantic Evaluation_ , pages 75–86, Minneapolis, Minnesota, USA. Association for Computational Linguistics.