\section{INTRODUCTION}
Object P/2010 H2 (Vales) was discovered by Jan Vales as a $m_V \sim$ 12.5 magnitude source in data taken with the 0.6 m telescope at Crni Vrh Observatory (Slovenia) on UT 2010 April 16.01 (Vales et al.~2010).  While first noticed as a point source, visual observers began within a day  to report an expanding coma (Vales et al.~2010, Mikuz et al.~2010).  Prediscovery observations include an upper  limit to the brightness  $m_V \gtrsim$ 20 on UT April 15.4 (Catalina Sky Survey observation by R. Kowalski, reported in Vales et al.~2010), another non-detection at unfiltered magnitude $>$16  on April 15.56 (Balanutsa et al.~2010) and a pre-discovery detection at $m_V$ = 13.7$\pm$0.1 on April 15.82 (Balanutsa et al.~2010), the latter two measurements with 0.4 m diameter telescopes of the MASTER network. Combined, the observations show sudden brightening of this previously unknown object with a doubling time of hours.  Outbursts in comets are not rare (e.g.~13 were tabulated in a 9 year period by Ishiguro et al.~2016), but few have amplitudes as large as the $\ge$ 7.5 magnitude range (a factor of $\ge 10^3$) in P/Vales, and few have been observationally well-characterized, providing a motivation for this study.

The orbit has semimajor axis, $a$ = 3.850 AU, eccentricity, $e$ = 0.193 and inclination, $i$ = 14.3\degr, an orbital period of 7.56 year and P/Vales passed perihelion (at 3.108 AU) on UT 2010 March 9.3.    The orbital elements give a Tisserand parameter with respect to Jupiter, $T_J$ = 2.99, signifying strong gravitational interactions with that planet.  This Tisserand is too small to qualify P/Vales as an active asteroid (Jewitt et al.~2015), for which $T_J \ge$ 3.08 is a minimum requirement.  P/Vales has instead been described as a  Hilda asteroid (Hildas are bodies in orbit near the location of the 3:2 mean-motion resonance with Jupiter, at semimajor axis $a$ = 3.971 AU).   P/Vales indeed lies within the  orbital element range of the Hilda group as defined  by Zellner et al.~(1985), namely 3.7 $\le a \le$ 4.2 AU, $e \le$ 0.3, and $i \le$ 20\degr.   However,  its  semimajor axis is smaller than 96\%, and its inclination larger than 92\%, of the $\sim$5000  objects in this semimajor axis range (Figure \ref{aei}).   Even though we refer to the osculating orbital elements (which change with time owing to perturbations from the planets), the figure serves to show that P/Vales is dynamically distinct from the typical Hilda asteroids clustered around the resonance.  Supporting this difference, Marsden (2010) found that P/Vales passed within 1 AU of Jupiter in 1976, and six comparably close encounters have occurred within the 20th century\footnote{JPL Small Body Database Browser: \url{https://tinyurl.com/y3kfyahe}}.  Consequently, we regard P/Vales as a likely ``Quasi-Hilda", presumably a temporarily captured Jupiter-family comet, several examples of which have been previously described in the literature (e.g.~Toth 2006, Gil-Hutton and Garcia-Migani 2016).

The outburst of P/Vales attracted considerable observational attention  over the spring and summer of 2010. We searched for additional detections of P/Vales  from the prior and subsequent orbits in archival data, finding none. While some of the  observations from 2010 were obtained by eye and are of limited scientific value, even  calibrated observations taken with electronic detectors  have failed to find their way into the refereed literature, a situation we begin to reverse with this paper.

\section{OBSERVATIONS}
Early epoch observations were taken at our request only 4 days after the discovery, on UT 2010 April 20, using the ``Baade'' 6.5 m diameter Magellan telescope, located in Chile, by Scott Sheppard. We used a 2048$\times$4096 pixel section of the IMACS short camera, giving a field of view approximately 400\arcsec$\times$800\arcsec~at the 0.2\arcsec~pixel$^{-1}$ image scale.  The seeing was 0.7\arcsec~to 0.8\arcsec~Full Width at Half Maximum (FWHM).    A total of 24 images of P/Vales were secured, with a range of exposure times from 2 s to 300 s.  After rejecting images in which the core of the target was saturated, we analysed 3 images in Sloan g', 9 in Sloan r', and 4 in Sloan i'.  Photometric calibration of the data was obtained using measurements of field stars, calibrated as part of the Deep Lensing Survey (DLS; Wittman et al.~2002, Smith et al.~2002).

Observations on UT 2010 June 6, 15, 18 and July 3 were obtained using the 0.9 m telescope of West Mountain Observatory, operated by Brigham Young University, in Utah.  We used the Finger Lakes PL-09000 charge-coupled device (CCD), which gives a 25.2\arcmin~field of view with 0.49\arcsec~pixels.  Images through a broadband Johnson-Cousins R filter were calibrated with reference to Landolt (1992) standard stars.

We also used the Keck 10 m diameter telescope atop Mauna Kea (altitude 4200 m) to observe P/Vales on UT 2010 August 20.  The Low Resolution Imaging Spectrometer (LRIS: Oke et al.~1995)  possesses independent blue and red channels separated by a dichroic filter.  We used the ``460'' dichroic (50\% peak transmission at 4900\AA~wavelength), and a broadband B filter on the blue side.  The B filter has central wavelength $\lambda_C$ = 4370\AA~ and Full Width at Half Maximum (FWHM) = 878\AA.  On the red side, we used a broadband R ($\lambda_C$ = 6417\AA, FWHM = 1185\AA) filter.   Photometric calibration of the data was secured using observations of  standard stars selected to have sun-like colors from the list by Landolt (1992) and cross-checked using field stars.  Unfortunately, the seeing was unusually poor and variable, in the range $\sim$2\arcsec~to 3\arcsec~FWHM.   As a result, we present only large aperture photometry from Keck, in order  to assess the total magnitude of P/Vales.

A timeline of observations is given in Table (\ref{geometry}) while representative composite images are shown in Figure (\ref{images}).

\subsection{Photometry}
We elected to measure photometry within a nested set of circular apertures having projected, fixed radii of 5$\times10^3$ km, 1$\times10^4$ km, 2$\times10^4$ km and 4$\times10^4$ km at the distance of P/Vales (Table \ref{photometry}).  The use of fixed linear (as opposed to angular) apertures ensures that we  measure the same volume around the nucleus independent of distance from the Earth, and obviates the need for a correction dependent on the surface brightness profile.  Sky subtraction was determined from the median signal measured in a contiguous annulus with inner radius 4$\times10^4$ km and extending out to (1 to 2) $\times 10^5$ km, depending on the dataset.    In some cases, we digitally removed background objects projected within the nested apertures.   We obtained photometric calibration of the data using   images of nearby Landolt (1992) standard stars and of field stars calibrated in the Pan STARRS survey (Tonry et al.~2012, Magnier et al.~2013).   Uncertainties on the photometry have several components.  The values listed in Table (\ref{photometry}) are the uncertainties determined empirically from the standard error on the mean of repeated measurements in each filter.  Small, additional uncertainties arise from the differences in the bandpasses  of the filters used at each telescope and, particularly for the Landolt stars, intrinsic uncertainties of the magitudes at the $\pm$0.01 to 0.02 magnitude level.  These  are all small compared to the uncertainty introduced by the unmeasured phase angle dependent darkening of P/Vales, and therefore of no consequence here.  All photometry was converted to standard Kron-Cousins BVR magnitudes using transformation equations from Smith et al.~(2002) for SDSS data and Tonry et al.~(2012) where Pan STARRS data were used.

We focus our analysis on the  apparent R-filter magnitudes, $m_{R}$, which we convert into absolute magnitudes, $m_R(1,1,0)$, using

\begin{equation}
m_{R}(1,1,0) = m_{R} - 5\log_{10}(r_H \Delta) - \beta_{ph} \alpha.
\label{H}
\end{equation}

\noindent Here, $r_H$ and $\Delta$ are the instantaneous heliocentric and geocentric distances expressed in AU, respectively, and $\alpha$ is the phase angle in degrees.  Quantity $\beta_{ph}$ is the phase coefficient, equal to the ratio of flux densities scattered at angle $\alpha$ to $\alpha$ = 0\degr.  The phase coefficient is unmeasured in P/Vales but studies of other comets show that in back-scattering (small $\alpha$) geometries,  $\beta_{ph}$ = 0.02 magnitudes degree$^{-1}$ provides a useful approximation (Meech and Jewitt 1987).

The absolute magnitude provides a measure of the scattering cross-section, $C$ (km$^2$), through

\begin{equation}
C = \frac{1.5\times10^6}{p} 10^{-0.4 m_{R}(1,1,0)}
\label{C}
\end{equation}

\noindent in which $p$ is the geometric albedo.  We assume $p$ = 0.1 throughout, consistent with the nominal albedo of cometary dust (Zubko et al.~2017) and with the range of albedos (0.03 to 0.12)  inferred in the possibly similar, outbursting comet 17P/Holmes (Ishiguro et al.~2010).  The surfaces of Hildas have lower average albedos, $p \sim$ 0.05 (Grav et al.~2012). As we note later, ice was reported in P/Vales (Yang and Sarid 2010), suggesting that a higher albedo might be appropriate.  On the other hand,  ice becomes optically absorbing and dark unless very pure.  In short,  it is not obvious what the albedo of the ejected material should be, or even that albedo should be constant with respect to time since ejection.  Cross-sections can be easily scaled to other albedos from Equation (\ref{C}) by the factor $0.1/p$.

\subsection{Color Photometry}

Color measurements  from UT 2010 April 20 are presented in Table (\ref{colors}), along with data for the colors of the Sun from Holmberg et al.~(2006).  All three independent color indices, B-V, V-R and R-I, show evidence for a trend towards smaller values at larger radii, indicating that the outer parts of the coma are more blue than the inner parts (Figure \ref{colorplot}).  It is unlikely that the gradients are caused by gas contamination of the signal, both because  the resonance fluorescence bands of common molecules are weak and rarely detected at 3 AU and because these bands are largely confined to  wavelengths $\lesssim$5000\AA~and would mainly affect B-V.  All three color gradients gradients are individually significant at the 4$\sigma$ level, and may be related to particle fragmentation inferred from the surface brightness profile (Section \ref{profile}).  However, we do not model this effect here, since color depends on many unknown properties of the dust (size, composition, microstructure), any or all of which may vary with time in a transient body like P/Vales, and unique models of optical colors cannot be constructed.  Instead, we use the colors for comparative purposes only.

The central colors are better representative of the source object and less likely to be affected by time-of-flight dependent optical effects in small particles.   Optically, the central (5000 km) colors of P/Vales are redder than the mean color of active short period comets (B-R = 1.40$\pm$0.02 vs.~1.22$\pm$0.02 for the comets; Jewitt 2015), but individual comets in the latter sample are scattered over the range B-R = 1.0 to 1.4, and some are as red as P/Vales.  The colors of P/Vales are much redder than the colors of  inactive Hildas. For example, the spectral slope across the B to R region is about $S'$ = 20\%/1000\AA, whereas the slopes measured for Hildas are all $S' \le$ 14\%/1000\AA~(Dahlgren and Lagerkvist 1995, Gil-Hutton and Brunini 2008). The red colors are more consistent with the mean colors of the nuclei of Jupiter family comets (for which B-V = 0.87$\pm$0.05, V-R = 0.50$\pm$0.03 and R-I = 0.46$\pm$0.03; Jewitt 2015).  This is not proof that P/Vales is a resonantly captured comet, but the colors are consistent with this interpretation.   A single large aperture (16.3\arcsec~radius, or 40,000 km) color measurement was obtained at Keck on UT 2010 August 10 under conditions of poor seeing, giving B-R = 1.29$\pm$0.10, consistent with the data from Table (\ref{colors}).

\subsection{Surface Brightness Profile}
\label{profile}
We measured the surface brightness profile, $\Sigma(\theta)$, where $\theta$ is the angular distance from the nucleus, in the UT 2010 April 20 data as follows.  The profile  was computed using a nested set of concentric annuli each 1 pixel (0.2\arcsec) wide, centered on the nucleus and with sky subtraction from a surrounding annulus having inner and outer radii 400 (80\arcsec) and 800 pixels (160\arcsec), respectively.  The profile was extracted out to a radius of 100 pixels (20\arcsec).  Similar measurements were taken to determine the point-spread function (PSF) from the profiles of field stars, using nearly simultaneous integrations of 10 s duration in which non-sidereal trailing is negligible.  The profiles are shown in Figure (\ref{sb}).

As expected, the coma profile shows a central excess caused by the convolution of the comet profile with the PSF.  At radius $\theta$ = 1\arcsec, however, the coma surface brightness is already $>10$ times larger than the PSF brightness and, at larger radii, the effects of the convolution can be ignored.   The profile becomes steeper at  angular radii $\theta \gtrsim$ 4\arcsec, as the physical edge of the coma is approached.  The logarithmic gradient of the surface brightness is $m = d\ln\Sigma(\theta)/d\theta$.  In the range 1\arcsec~$\le \theta \le$ 4\arcsec~(1545 to 6180 km at the comet), we find $m = -0.83\pm$0.01.  The gradient is significantly less steep than the canonical $m$ = -1 expected of a steady-state coma from the equation of continuity (Jewitt and Meech 1987).    This could be because the coma is itself not in steady state, with production at the nucleus varying on a timescale comparable to or shorter than the residence time for dust particles.  Alternatively,  the shallow gradient could reflect fragmentation of ejected particles, resulting in the creation of extra scattering cross-section (and perhaps a change in broadband color) as distance from the nucleus increases.  This was the  case in 17P/Holmes, where an even flatter surface brightness profile, $m$ = -0.27, was recorded (Stevenson and Jewitt 2012).  The fractional increase in the cross-section varies as $\theta^{1+m} - 1$.  With $\theta$ = 4\arcsec~and $m$ = -0.83, for example, the increase in the cross-section required to fit the gradient is  by a modest  $\sim$27\%.

The spectroscopic detection of water ice in the comae of both 17P/Holmes (Yang et al.~2009) and P/Vales (Yang and  Sarid 2010)  suggests a role for sublimation.  Specifically,  composite grains bound together by water ice would spontaneously disaggregate  through sublimation following their sudden expulsion from the nucleus into sunlight.  The ice sublimation rate is a strong function of the grain temperature as set by the heliocentric distance and albedo.  We solved the energy balance equation for an exposed, sublimating ice surface (c.f.~section 4.1 of Jewitt et al.~2020).  As an example, at $r_H$ = 3.11 AU, we calculate that the sublimation rate from an isothermal water ice sphere of albedo 0.1 is $f_s$ = 2$\times10^{-8}$ kg m$^{-2}$ s$^{-1}$.  The resulting sublimation lifetime of a grain of radius $a$ is $t_s \sim a \rho/f_s$.  By substitution, we estimate $t_s \sim$ 0.6 to 6 day, for $a$ = 1 to 10 $\mu$m.  Since $t_s$ is  comparable to the $\sim$4 day interval between the start of the outburst and the April 20 observation, we consider sublimation-induced disaggregation to be a plausible explanation of the shallow surface brightness gradient.  Spatially and/or temporally resolved spectroscopic observations of ice absorption could test this explanation.

\section{DISCUSSION}
\subsection{Analytic Considerations}
\label{simple}

\textit{Nucleus Radius, $r_n$:} The strongest observational constraint on the nucleus radius is set by the non-detection at $m_V \ge$ 20 from the Catalina Sky Survey on UT 2010 April 15.4, when $r_H$ = 3.112, $\Delta$ = 2.131, $\alpha$ = 4.7\degr.   By Equation (\ref{H}) these values give absolute magnitude $m_V(1,1,0) \ge$ 15.8 and, by Equation (\ref{C}), a nucleus cross-section $C_n \le$ 7.2 km$^2$ and effective radius $r_n = (C_n/\pi)^{1/2} \le$ 1.5 km.

\textit{Start of Outburst, $T_0$:}  The non-detections on UT 2010 April 15.4 (Vales et al.~2010) and 15.56 (Balanutsa et al.~2010) and the first detection on April 15.82 strongly bracket the beginning of the outburst.  To obtain a better estimate of the start time within this range, we summarize the early-time observations in Figure (\ref{rise}), where we have converted the reported apparent magnitudes and limits to cross-sections, as described above.   Upper limits to the cross-sections are marked in the figure with down-pointing arrows.  We fitted the three detections with an exponential function $C = C_0(1- s \exp(-t/w))$, where $C_0$ = 11,069 km$^2$ is the peak cross-section determined on UT 2010 April 20.2 and $s$ and $w$  are constants.   The function, plotted in the figure, extrapolates to $C$ = 0 on UT 2010 April 15.76 (DOY 105.76) and provides our best estimate of the time of initiation.  This is only 0.2 days (5 hour) after the non-detection reported by Balanutsa et al.~(2010), 0.06 days (1.4 hour) before the first detection by the same observers and 0.25 days (6 hour) before the discovery.  Thus, even though the adoption of an exponential function is arbitrary, the data are highly constraining and there is little room for the start-time to be much different from that indicated by the fit.  The e-folding rise time of the cross-section given by the fit is $w^{-1}$ = 0.22 days (5.4 hour), corresponding to a rise-time half-life $t_{1/2} = \ln(2)w^{-1}$ = 3.7 hour.  Even shorter rise-times (0.3 to 0.6 hour) were reported in early-phase observations of 17P/Holmes (Hsieh et al.~2010).  The peak measured cross-section, $C_{max} = 1.8\times10^4$ km$^2$ on UT 2010 April 20.2 (Table \ref{photometry}), is equal to that of a circle of radius  $r \sim (C_{max}/\pi)^{1/2}$, or $r \sim$ 75 km.

\textit{Particle Parameters:} The detailed appearance of the coma is influenced by the time-profile of the emission and by the distribution of particle radii. We take the latter to be a differential power law with index, $q$, such that the number of particles with radius between $a$ and $a+da$ is $n(a) da = \Gamma a^{-q} da$ in the range $a_{min} \le a \le a_{max}$, with $\Gamma$ equal to a constant.  Before we apply a Monte-Carlo approach to constrain the dust parameters, it is informative to use basic physics to assess the nature of the particles.  We next obtain values for $a_{min}$, $a_{max}$ and $q$ to compare with values determined independently from a numerical model.

\textit{Lower Size Limit $a_{min}$:} We make two morphological observations from the images taken on UT 2010 April 20.2.   First, the coma extends in the sunward direction by a distance $\ell \sim 4\times10^7$ m.  We interpret $\ell$ as the turn-around distance for dust particles ejected sunward at speed $U$, and subject to a constant radiation pressure induced acceleration.  We write this acceleration as $\beta g_{\odot}$, where radiation pressure efficiency factor $\beta$ is a dimensionless number and $g_{\odot} = GM_{\odot}/r_H^2$ is the local gravitational acceleration to the Sun.   Parameter $\beta$ is a function of particle size, shape and composition.  As a useful approximation for dielectric spheres, we take $\beta \sim 10^{-6}/a$, where $a$ is the particle radius expressed in meters (Bohren and Huffman 1983).  Then, the equation of motion for a fixed acceleration gives

\begin{equation}
\beta  = \frac{2\ell r_H^2}{G M_{\odot} t^2}.
\label{beta}
\end{equation}

\noindent where  $G = 6.67\times10^{-11}$ N kg$^{-2}$ m$^2$ is the gravitational constant, $M_{\odot} = 2\times10^{30}$ kg  is the mass of the Sun and $t$ is the time of flight.

Second, the comet possesses no clear radiation-swept tail (i.e.~dust particles launched sunward have been propelled to the anti-solar side of the nucleus by a distance no greater than $-\ell$).   Again, from the equation of motion for a fixed acceleration of $\beta g_{\odot}$, and substituting for $\beta$ from Equation (\ref{beta}),  we find

\begin{equation}
U = \frac{2\ell}{t}
\label{U}
\end{equation}

\noindent For the Magellan observation we set $t$ = 4.45 days (3.8$\times10^5$ s) and $\ell = 4\times10^7$ m to find $U$ = 210 m s$^{-1}$ and $\beta$ = 0.9 (implying $a \sim$ 1.1 $\mu$m).  Strictly, $\ell$ is a lower limit to the true turn-around distance because of the effects of projection into the plane of the sky. Therefore, the derived values of $U$ and $\beta$ are also lower limits.  For present purposes, however, these crude estimates are sufficient to show that  the early-stage morphology of the envelope of the coma  is controlled by fast-moving particles with radii $a \sim$ 1 $\mu$m.

By UT 2010 June 6 ($t$ = 52 days after the outburst), the tail of P/Vales can be traced to the edge of the field of view, a distance $\ell \ge 1.3\times10^9$ m (Figure \ref{jun06}).  Application of Equation (\ref{beta}) then gives $\beta \ge$ 0.3 for the dust particles at the edge of the field of view.  Again, the large-scale morphology on this date is controlled by  the dynamics of micron-sized particles accelerated by  solar radiation pressure.  Large particles, less easily accelerated by radiation pressure,  must be located closer to the nucleus on this date.  The absence of a distinct feature attributable to large particles (specifically, a dust trail in data from UT 2010 June 6 (Figure \ref{images}), when the Earth was only 1.1\degr~above the orbital plane of P/Vales (Table \ref{geometry})), shows that the large particle contribution to the cross-section of the ejected material is minor on this date.

\textit{Upper Size Limit $a_{max}$:} Larger particles are less well-coupled to the outrushing gas and will have smaller terminal speeds.  In gas drag expulsion, the terminal speed of a dust grain is related to its radius by $U(a) = U_1(a_0/a)^{1/2}$, where constant $U_1$ is the speed of a particle at reference radius $a_0$.  We take $a_0$ = 10$^{-6}$ m and $U_1$ = 210 m s$^{-1}$.  To estimate the radius of the largest (slowest) particle that can be ejected against the gravity of the nucleus (ignoring possible effects of cohesion), we set  $U = V_e$, where $V_e$ is the gravitational escape speed, to find

\begin{equation}
a_{max} = \frac{3 a_0 U_1^2}{8 \pi G \rho r_n^2}
\label{amax}
\end{equation}

\noindent where  we have assumed for simplicity that the nucleus is spherical and of density $\rho$. Substituting $r_n $ = 1.5 km, we find from Equation (\ref{amax})  the largest ejectable dust radius,  $a_{max} \sim$ 6 cm.  These large, low speed particles would take a time $t \sim  3\times10^7$ s (1 year) to cross the 40,000 km photometry aperture, consistent with the persistence of excess cross-section in the lightcurve (Figure \ref{Ce_vs_DOY}) on 100 day timescales.

\textit{Size Distribution Index, $q$:} Smaller, faster dust particles escape the photometry aperture more quickly than larger, slower ones.  Therefore, in an impulsive outburst, the mean size of the particles in the aperture should increase with time because of the preferential loss of small particles, even as the total cross-section decreases.  An expression for the fading caused by the selective loss of particles from a photometry aperture was derived by  Jewitt et al.~(2017).  Here, we modify this expression to consider particles within an annulus, rather than a circular aperture, in order to compare with annular data extracted from Table (\ref{photometry}).

Consider a photometry annulus with inner and outer radii, $r_1$ and $r_2$, respectively, and a measurement taken time $t$ after the ejection.  Particles traveling more slowly than $U_{min} = r_1/t$ will not have reached the inner edge of the annulus, while those traveling faster than $U_{max} = r_2/t$ will have escaped its outer edge.   Substituting for the velocity-size relation, $U(a) = U_1(a_0/a)^{1/2}$,  we find the critical radii

\begin{equation}
a_1(t) = a_0 \left(\frac{U_1 t}{r_1}\right)^2~~\textrm{and}~~a_2(t) = a_0 \left(\frac{U_1 t}{r_2}\right)^2.
\end{equation}

\noindent Particles with radii $a < a_2(t)$ will have escaped the outer edge of the annulus after time of flight, $t$, while those with $a > a_1(t)$ travel so slowly that they have not yet reached the inner edge.  Therefore, at any time, $t$, the particles within the annulus are  confined to the radius range $a_2(t) < a < a_1(t)$, assuming impulsive emission

The fraction of the ejected dust cross-section remaining within the annulus at time $t$ since ejection is then

\begin{equation}
\frac{\Delta C(t)}{C(0)} = \frac{\int_{a_2}^{a_1}\pi a^2 n(a) da}{\int_{a_{min}}^{a_{max}}\pi a^2 n(a) da}
\end{equation}

\noindent where $C(0)$ is the total ejected cross-section.  Substituting $n(a)da = \Gamma a^{-q}da$ and evaluating, we obtain

\begin{equation}
\frac{\Delta C(t)}{C(0)}  = \left(\frac{a_0}{a_{min}}\right)^{3-q}  \left[\left(\frac{1}{r_1}\right)^{6-2q} - \left(\frac{1}{r_2}\right)^{6-2q}\right]   (U_1 t)^{6-2q}
\label{model}
\end{equation}

\noindent where we have assumed $a_{max} \gg a_{min}$, $q > 3$ and that the size range $a_1 - a_2$ is contained with $a_{min} - a_{max}$.  The  time dependence in Equation (\ref{model}), $\Delta C(t)/C(0) \propto t^{6-2q}$, is the same as derived in Jewitt et al.~(2017), and allows us to estimate $q$ from the time-dependent annular photometry.

The four annular cross-sections, $\Delta C(t)$, are plotted in Figure (\ref{comafit}).  We least-squares fitted power laws to $\Delta C(t)$  in order to find the index, $q$ from Equation (\ref{model}).  Lines in the figure show independent fits to photometry from the four apertures presented in Table (\ref{photometry}).  For the $\Delta C$ = 0-5,000 km, 5,000-10,000 km, 10,000-20,000 km, 20,000-40,000 km apertures we find, respectively, $q$ = 3.80, 3.70, 3.61 and 3.50.  The formal uncertainties on these fits are $\pm$0.01 to $\pm$0.02, but we use the standard error on the mean of the four values as a better measure of the true scatter.  The mean size distribution index from fits to all four apertures is $q$ = 3.61$\pm$0.06.  Small differences between the values could result from many causes (perhaps the size distribution is not a power law, as assumed in Equation (\ref{model}), perhaps the emission was not impulsive or the size-speed relation is invalid, perhaps the grain albedo changes with time as ice sublimates away).  Given these many potential problems, the overall  agreement between the four independent determinations of $q$ is good.  The derived value is close to $q$ = 3.7$\pm$0.1 determined for fragmenting asteroid P/2012 F5 (Gibbs) (Moreno et al.~2012) and $q$ = 3.6$\pm$0.6 for fragmenting comet 332P/Ikeya-Murakami (Jewitt et al.~2016).

\textit{Dust Mass, $M$:} The mass of an optically thin collection of spheres is related to the sum of their cross-sections by

\begin{equation}
M = \frac{4}{3}\rho \overline{a} C
\label{mass}
\end{equation}

\noindent  where $\overline{a}$ is the area-weighted mean particle size responsible for cross-section, $C$.  The latter is given by

\begin{equation}
\overline{a} = \frac{\int_{a_{min}}^{a_{max}}a^3 \Gamma a^{-q} da}{\int_{a_{min}}^{a_{max}}a^2 \Gamma a^{-q} da}
\end{equation}

\noindent which, for $q$ = 3.61 simplifies to

\begin{equation}
\overline{a} = 1.56 a_{max}^{0.39} a_{min}^{0.61}
\label{abar}
\end{equation}.

\noindent Substitution gives $\overline{a} = 1\times10^{-4}$ m (100 $\mu$m).

We see from  Table (\ref{photometry}) and Equations (\ref{mass}) and (\ref{abar})  that the peak dust cross-section, $C$ = 17.6$\times10^3$ km$^2$, with $\rho$ = 500 kg m$^{-3}$, gives dust mass $M = 1.2\times10^{9}$ kg.   This is about 10$^{-4}$ times the mass of the nucleus, considered as a sphere of the same density and radius 1.5 km.  The bulk of this mass was released over a period $\tau \le$ 1 day (Figure \ref{rise}), corresponding to a  mean mass production rate $M/\tau \gtrsim$ 13,600 kg s$^{-1}$.   This rate is comparable to the (1 to 2)$\times10^4$ kg s$^{-1}$ measured in super-active comet C/1995 O1 (Hale-Bopp) at similar heliocentric distance (Weiler et al.~2003), although it is sustained in P/Vales for only a short period of time.

\textit{Specific Energy, E/M:}
The energy per unit mass of ejecta is $E/M = (1/2) U^2(\overline{a})$.  For a $q$ = 3.6 distribution, the energy is carried by the small (fast) particles, because of the $U^2$ dependence, while the mass is carried by the large particles.  With $\overline{a}$ = 100 $\mu$m, $U(\overline{a}$) = 21 m s$^{-1}$ and $E/M$ = 220 J kg$^{-1}$.  This specific energy is about 10$^2$ times smaller than the corresponding quantity estimated for comets 332P/Ikeya-Murakami and 17P/Holmes.  The total energy of the outburst is $E \sim 3\times10^{11}$ J, an order of magnitude smaller than the outburst of 332P/Ikeya-Murakami and at least 10$^3$ times smaller than that of 17P/Holmes (summarized by Ishiguro et al.~2016).

\subsection{Monte-Carlo Model}
\label{MC}
We applied a 3D Monte-Carlo simulation
(Ishiguro et al.~2007) to P/Vales. In the model the motions of dust particles,  under the action of solar gravity and radiation pressure, are followed as a function of particle size and time and direction of ejection.   Our aim is not to reproduce the comet exactly since, given the number of free parameters in the model, a match can almost always be obtained.   Instead, we aim to provide a consistency check of the results deduced analytically in Section (\ref{simple}).

As in Section (\ref{simple}), we set the dust velocity-size relation, $U(a) = U_1(a_0/a)^{1/2}$, as expected of gas drag acceleration, with $U_1$ = 210 m s$^{-1}$.  We further assumed the range of particle sizes 10$^{-6} \le a \le 6\times10^{-2}$ m, a differential power-law index $q$ = 3.6 and impulsive ejection on UT 2010 April 15.76.    The model includes, as a free parameter, the angular dependence of the dust production rate, represented by  a cone having its apex at the center of the nucleus, the Sun on its axis and a variable half-width, $w$. By experimentation, we found that   $w$ is most strongly constrained by the angular width of the dust tail in data from 2010 June, and that there is a trade-off between $w$ and $U_1$ in the dust speed vs.~size relation.  Given the nominal speed $U_1$ = 210 m s$^{-1}$, we found that $w$ = 25\degr~provided an acceptable fit to the tail.  Values $w$ = 15\degr~and $w$ = 35\degr~produced tails, respectively, too narrow and too wide to fit the images from 2010 June.

Examples of the Monte Carlo simulations are shown in Figure (\ref{mc}), where it is seen that the parameters deduced analytically indeed generate models that match the data well.
For example, the  simulations  show that, with the parameters deduced in Section (\ref{simple}), no large-particle trail can be discerned even near the crossing of the orbit plane on UT 2010 June 6 (Table \ref{geometry}).  This is consistent with the absence of a trail in the data (Figure \ref{mc}), confirming that large, slowly ejected grains present only a small fraction of the total cross-section even when spatially dispersed in the plane of the sky by radiation pressure. The simulations also reproduce the persistence of near-nucleus dust into 2010 August and the roughly circular late-stage appearance of the comet (Figure \ref{mc}).  This confirms that, by August 10,  only the largest, slowest particles remain in the vicinity of the nucleus and the smaller particles, constituting the elongated tail in data from 2010 June, have been swept away.  Overall, Monte Carlo simulations based on  the parameters deduced in Section \ref{simple} provide an acceptable match to the broad-brush appearance of P/Vales.

Closer examination of the data shows evidence for unmodeled anisotropy in the coma in addition to that caused by solar radiation pressure.  In Figure (\ref{oneoverrho}) we show, on the left, the image from UT 2010 April 20 and, on the right, the same image divided by a normalized profile in which the surface brightness varies inversely with angular distance from the nucleus (the location of which is marked in the figure by a black dot; Samarasinha and Larson 2014).  The filtered image clearly shows excess emission near position angle 220\degr, and a broad deficit towards 330\degr.  The amplitudes of these features are small (explaining why they are not evident in the left, unfiltered, image) but they confirm that mass loss from P/Vales was not isotropic.  We have not attempted to model the anisotropy evident in  Figure (\ref{oneoverrho}).

\subsection{Outburst Mechanism}  We briefly consider possible causes of the outburst of P/Vales.

The existence of high speed ejecta, with $U \ge 210$ m s$^{-1}$ for 1 $\mu$m particles,  strongly suggests the action of gas drag.  To see this, we note that, at $r_H$ = 3.112 AU, the local isothermal blackbody temperature is $T_{BB}$ = 158 K while the mean thermal speed of H$_2$O molecules at this temperature is $V_{th} = (8 k T_{BB}/(\pi \mu m_H))^{1/2}$, where $\mu$ = 18 is the molecular weight of water and $m_H = 1.67\times10^{-27}$ kg is the mass of the hydrogen atom.  We compute $V_{th}$ = 430 m s$^{-1}$.  The observation that, within a factor of two, $U \sim V_{th}$ for micron-sized dust grains is consistent with their acceleration by gas drag, as is the case in comets, generally.

\textbf{Volatile Explosion:}
Evidence for the interplay between cometary activity and dynamics is well-established.  For example, the largest values of the non-gravitational accelerations of comets (which provide a measure of the mass loss rate) are correlated with recent,  inward migration of the perihelion distance (Rickman et al.~1991).   Could it be that the recent trapping of P/Vales into its quasi-Hilda orbit has triggered either enhanced sublimation of sub-surface volatiles or the exothermic crystallization of buried amorphous ice, leading to pressure build-up and eventual rupture of  a cohesive mantle (c.f.~Samarasinha 2001)?

In this scenario, the  delay between the  Jupiter encounters of the previous century and the  outburst in 2010 would  be a consequence of the slow conduction of heat from the surface to the buried ice.   Solution of the conduction equation shows that the distance, $d$, over which heat is conducted in time $\tau$ is approximately $d \sim (\kappa \tau)^{1/2}$, where $\kappa$ is the thermal diffusivity of the material.  The diffusivity of solid dielectric materials is $\kappa \sim 10^{-6}$ m$^2$ s$^{-1}$, but porosity reduces $\kappa$ significantly. For example, laboratory measurements of highly porous dielectric powders give $\kappa \sim 10^{-9}$ to $10^{-8}$ m$^2$ s$^{-1}$ (Sakatani et al.~2018).  If we suppose that the most recent inward excursion of the perihelion distance occurred in 1973, the resulting $\tau =$ 35 years (10$^9$ s) delay would correspond to a buried ice depth $d \sim$ 1 m to 3 m.  Earlier perihelion decreases would correspond to larger $\tau$ and deeper ice, but the diffusive $d \propto \tau^{1/2}$ dependence suggests that  $d$ is unlikely to be substantially larger than $\sim$10 m.

The Monte-Carlo models (Section \ref{MC}) indicate ejection into a cone of half-angle $w \sim$ 25\degr, amounting to a solid angle $\Omega  \sim$ 0.2 steradian.   Given this, we estimate the area of the region producing the outburst as $A_{ice} = 4\pi r_n^2 (\Omega /4\pi)$  and the thickness of ice needed to provide the ejecta mass, $M$, as $d_{ice} \sim M/(\rho \Omega r_n^2)$.  Setting  $M = 1.2\times10^9$ kg, $\rho$ = 500 kg m$^{-3}$ and $r_n \le 1.5$ km, we find  $A_{ice} \le$ 0.4 km$^2$ (1.5\% of the nucleus surface) and $d_{ice} \ge$ 5 m.  Thus, a volatile pocket  modest in both thickness and areal extent  could supply the ejecta responsible for the outburst in P/Vales.    The crystallization of amorphous water ice releases $\sim10^5$ J kg$^{-1}$, far more than the average $\sim$220 J kg$^{-1}$ measured for the specific energy of the ejecta.   Thus, even a very inefficient conversion of crystallization energy into kinetic energy would be sufficient to drive the outburst and we consider a delayed explosion caused by slowly conducted heat to be a plausible explanation for the outburst of P/Vales.  On the other hand, whether or not amorphous ice could persist at meter depths in P/Vales depends on its past dynamical history which is, as with all comets,  chaotic and unknown.    We note that the spectrum reported by Yang and Sarid (2010) showed the 1.65 $\mu$m wavelength absorption indicative of  water ice in the crystalline state.  However, at the local blackbody temperature (158 K at 3.1 AU), crystallization occurs in $\ll$1 second, and the initial state of the ice before ejection into sunlight cannot be spectroscopically ascertained.

\textbf{Cliff Collapse:}   Could the outburst of P/Vales have resulted from the collapse of a cliff or an overhang on the nucleus surface?   In this scenario, the collapsed material crumbles and spreads across the surface as a landslide from which the entrained volatiles  sublimate, expelling dust and debris.  A high-temperature Terrestrial analog might be found in the pyroclastic flows produced by gas-rich, high viscosity volcanic magmas, where ejected solids degas and move as a ground-hugging, dense flow.  A cliff-collapse outburst of $\sim10^6$ kg  has been observed close-up on the nucleus of 67P/Churyumov-Gerasimenko (Pajola et al.~2017).

While qualitatively attractive, cliff collapse offers a less obvious explanation of the much larger ($\sim10^9$ kg) outburst in P/Vales.  To see this, we first note that the area $A_{ice}$, estimated above, corresponds in this scenario to the area of the landslide needed to generate the outburst by sublimation.
We solved the energy balance equation for sublimating carbon monoxide (CO) ice, considered as a representative supervolatile, finding that at $r_H$ = 3.112 AU, the specific sublimation rate of a perfectly absorbing CO ice surface oriented normal to the Sun is $f_s(CO) = 4.6\times10^{-4}$ kg m$^{-2}$ s$^{-1}$.   The sublimation rate from the landslide is then just $dM/dt = f_s(CO) A_{ice}$.  By substitution, $dM/dt \sim 180$ kg s$^{-1}$.  This is an upper limit to the true $dM/dt$ because the ice is unlikely to be perfectly absorbing (Pajola et al.~2017 determined an albedo $\sim$0.4 in their example) and because the landslide surface  is unlikely to be oriented perpendicular to the direction of the Sun.  Significantly, $dM/dt$ from this estimate is two orders of magnitude smaller than the  actual production rate obtained from the rising lightcurve,  $M/\tau = 13,600$ kg s$^{-1}$ (Section \ref{simple}). Sublimation from a fresh landslide produced by cliff collapse cannot directly supply the massive outburst of P/Vales.

However, expulsion of ice-containing debris from the landslide surface into an expanding, avalanche-like cloud having cross-section $\gg A_{ice}$ could strongly amplify the production rate.  For example, we found that the peak measured cross-section is equal to that of a circle of radius  $r \sim (C_{max}/\pi)^{1/2}$, or $r \sim$75 km (Section \ref{simple}).     Ice in these particles would then sublimate at a rate $\sim \pi r^2 f_s(CO)$ which, with $r$  = 75 km, could supply $dM/dt \sim 8\times10^6$ kg s$^{-1}$, two orders of magnitude larger than required by the data.   Furthermore, if the cliff collapse were to expose amorphous ice, the resulting immediate exothermic crystallization could easily supply enough energy to drive the outburst.  Particles with the mean radius $\overline{a}$ = 100 $\mu$m have velocity $U \sim$ 20 m s$^{-1}$, and would take $\sim$1 hour to fill an optically thin hemisphere of radius 75 km, a timescale consistent with the short rise-time of the lightcurve (Figure \ref{rise}).  The coma was probably optically thin even at the time of the first detection (April 15.82, about 0.06 days after the best-fit initiation time on April 15.76, corresponding to a delay of $\sim$1.5 hours).  We conclude that cliff collapse of a sufficient volume could generate an outburst having the magnitude and the rapid rise-time implied by the observations.

An unresolved issue with the cliff collapse hypothesis is one of timing.  Why, other than by coincidence, would cliff collapse occur decades after the entrapment of P/Vales into the 3:2 mean-motion resonance?  Perhaps low-level activity driven by the sublimation of near-surface water ice in the years prior to outburst caused incipient instability leading to collapse.  In the absence of relevant observational evidence, however, we can say nothing about this possibility.

\textbf{Impact:}
High speed ejecta can also be produced by impact.  However, two factors make an impact origin unlikely.  First, the collision probability in the Hilda population ($P_i \sim 2\times10^{-18}$ km$^{-2}$ yr$^{-1}$; Dahlgren 1998, dell'Oro et al.~2001)  is lower than in the main-asteroid belt ($P_i \sim 5\times10^{-18}$ km$^{-2}$ yr$^{-1}$, Bottke and Greenberg 1993). The Hilda population is $\lesssim$1\% of  the main-belt population and the quasi-Hilda population is orders of magnitude smaller still.  Second, and more seriously, the quasi-Hildas have short dynamical lifetimes in the 3:2 resonance region (e.g.~10$^3$ - 10$^4$ year, Gil-Hutton \& Garc{\'\i}a-Migani 2016). The likelihood of a substantial  impact in this small population with a tiny residence window is very small.  In the particular case of P/Vales, impact within a few decades, or even centuries, of its  injection into the present orbit (Marsden 2010) is incredibly unlikely.

\textbf{Other Processes: } The active asteroids are driven by several low energy processes in addition to sublimation gas-drag and impact (Jewitt 2012, Jewitt et al.~2015).  These include rotational disruption of the parent body, cracking caused by thermal expansion and/or desiccation stresses, and electrostatic ejection of fine dust.  None of these processes can offer a convincing explanation for  the high speeds reached by the ejecta in P/Vales.   In rotational disruption, for example, released material escapes with approximately the equatorial velocity of the rotation, typically $\le$ 1 m s$^{-1}$ at breakup for a kilometer-sized body.  Thermal fracture,  with realistic efficiencies for conversion of strain energy into kinetic energy of ejected fragments in the 1\% to 10\% range, produces peak speeds  $\sim$ (1 to 5) m s$^{-1}$ (Equation 14 of Jewitt 2012).  Electrostatic forces are weak, eject particles at similarly low speeds, and  cannot launch particles greater than micron-sized.  The mismatch with the ejecta speeds measured in P/Vales effectively  eliminates rotational disruption, thermal stresses and electrostatics as relevant mechanisms.

\clearpage

\section{SUMMARY}
We present time-resolved  observations of the 2010 photometric outburst of  P/2010 H2 (Vales). This is a quasi-Hilda  object, probably emplaced near the 3:2 mean-motion resonance  following a series of close encounters with Jupiter in the previous century.

\begin{enumerate}

\item The outburst, by $\Delta m_V\ge$ 7.5 magnitudes, started within a few hours of UT 2010 April 15.76, when 37 days past perihelion and at heliocentric distance 3.112 AU.

\item Ejected particles had a maximum cross-section 17,600 km$^2$ and mass $\sim 1.2\times10^9$ kg (10$^{-4}$ of the nucleus mass), with radii from microns to centimeters following a differential size distribution $n(a)da \propto a^{-3.61\pm0.06}$. Peak ejection rates were 13,600 kg s$^{-1}$, and  the ejection appears to have been impulsive, with a timescale $\lesssim$1 day.

\item High measured   particle ejection speeds (up to 210 m s$^{-1}$) are compatible with gas drag acceleration from sublimated ice.
They are incompatible with  rotational instability, thermal and desiccation stress fracture, and electrostatic repulsion, all of which are therefore ruled out as driving mechanisms.   Asteroid impact, while  capable of generating high-speed ejecta, is improbable given the small population and short dynamical lifetime of the quasi-Hildas.

\item P/Vales is most likely a temporarily captured comet in which conductive heating of   sub-surface ice has triggered an outburst, perhaps through  exothermic crystallization from the  amorphous state.  If so, an ice volume about 0.4 km$^2$ in areal extent (about 1\% of the nucleus surface) and $\ge$5 m thick, and buried beneath a refractory layer a few meters thick, is inferred.
\end{enumerate}

\acknowledgments
We thank Scott Sheppard for data from Magellan and comments on the manuscript, and Eric Hintz and Michael Joner  for data from the West Mountain Observatory.  We also thank Bin Yang and the two anonymous referees for comments. Some of the data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and NASA. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation.  Y.K. was supported by the European Research Council (ERC) Starting Grant No. 757390 (CAstRA).

{\it Facilities:}  \facility{Magellan Observatory, Keck Observatory, West Mountain Observatory}.


\section{Introduction}
\label{sec:intro}
Deeper and larger neural networks have led to significant improvement in accuracy in various tasks, but they are also more computationally expensive and unfit for resource-constrained scenarios such as online serving. An interesting and viable solution to this problem is knowledge distillation (KD) \citep{Bucilua:2006:MC:1150402.1150464,NIPS2014_5484,44873}, which can be used to transfer the knowledge of a large model (the teacher) to a smaller model (the student).
In the field of natural language processing (NLP), for example, KD has been successfully applied to compress massive pretrained language models such as BERT \citep{devlin-etal-2019-bert} and XLM-R \citep{conneau-etal-2020-unsupervised} into much smaller and faster models without significant loss in accuracy \citep{tang2019distilling,sanh2019distilbert,tsai-etal-2019-small,mukherjee-hassan-awadallah-2020-xtremedistil}.

A typical approach to KD is letting the student mimic the teacher model's output probability distributions on the training data by using the cross-entropy objective.
For structured prediction problems, however, the output space is exponentially large, making the cross-entropy objective intractable to compute and optimize directly.
Take sequence labeling for example. If the size of the label set is $L$, then there are $L^n$ possible label sequences for a sentence of $n$ words and it is infeasible to compute the cross-entropy by enumerating the label sequences.
Previous approaches to structural KD either choose to perform KD on local decisions or substructures instead of on the full output structure, or resort to Top-K approximation of the objective \citep{kim-rush-2016-sequence,kuncoro-etal-2016-distilling,wang-etal-2020-structure}.

In this paper, we derive a factorized form of the structural KD objective based on the fact that almost all the structured prediction models factorize the scoring function of the output structure into scores of substructures. If the student's substructure space is polynomial in size and the teacher's marginal distributions over these substructures can be tractably estimated, then we can tractably compute and optimize the factorized form of the structural KD objective. As will be shown in the paper, many widely used structured prediction models satisfy the assumptions and hence are amenable to tractable KD. In particular, we show the feasibility and empirical effectiveness of structural KD with different combinations of teacher and student models, including those with incompatible factorization forms.
We apply this technique to structural KD between sequence labeling and dependency parsing models under four different scenarios

\begin{enumerate} [leftmargin=*]
    \item The teacher and student share the same factorization form of the output structure scoring function.
    \item The student factorization produces more fine-grained substructures than the teacher factorization.
    \item The teacher factorization produces more fine-grained substructures than the student factorization.
    \item The factorization forms from the teacher and the student are incompatible.
\end{enumerate}

In all the cases, we empirically show that our structural KD approaches can improve the student models. In the few cases where previous KD approaches are applicable, we show our approaches outperform these previous approaches.
With unlabeled data, our approaches can further improve student models' performance.
In a zero-shot cross-lingual transfer case, we show that with sufficient unlabeled data, student models trained by our approaches can even outperform the teacher models.

\section{Background}

\subsection{Structured Prediction}
\label{sec:str_pred}
Structured prediction aims to predict a structured output such as a sequence, a tree or a graph. In this paper, we focus on structured prediction problems with a discrete output space, which include most of the structured prediction tasks in NLP (e.g., chunking, named entity recognition, and dependency parsing) and many structured prediction tasks in computer vision (e.g., image segmentation). We further assume that the scoring function of the output structure can be factorized into scores of a polynomial number of substructures. Consequently, we can calculate the conditional probability of the output structure ${\bm{y}}$ given an input ${\bm{x}}$ as follows:
\begin{align}
P({\bm{y}}|{\bm{x}}) & = \frac{\exp{(\text{Score}({\bm{y}},{\bm{x}})})}{\sum_{{\bm{y}}^{\prime} \in {\mathbb{Y}}({\bm{x}})} \exp{(\text{Score}({\bm{y}}^{\prime},{\bm{x}})})}\nonumber \\
            & = \frac{\prod_{{\bm{u}} \in {\bm{y}}}\exp{(\text{Score}({\bm{u}},{\bm{x}})})}{\mathcal{Z}({\bm{x}})}
\label{eq:2.1}
\end{align}
where ${\mathbb{Y}}({\bm{x}})$ represents all possible output structures given the input ${\bm{x}}$, $\text{Score}({\bm{y}},{\bm{x}})$ is the scoring function that evaluates the quality of the output ${\bm{y}}$, $\mathcal{Z}({\bm{x}})$ is the partition function,
and ${\bm{u}} \in {\bm{y}}$ denotes that ${\bm{u}}$ is a substructure of ${\bm{y}}$. We define the substructure space
${\mathbb{U}}({\bm{x}}) = \bigcup_{{\bm{y}}\in{\mathbb{Y}}({\bm{x}})}\{{\bm{u}}|{\bm{u}}\in{\bm{y}} \} $as the set of substructures of  all possible output structures given input ${\bm{x}}$.

Take sequence labeling for example. Given a sentence ${\bm{x}}$ , the output space ${\mathbb{Y}}({\bm{x}})$ contains all possible label sequences of ${\bm{x}}$. In linear-chain CRF, a popular model for sequence labeling, the scoring function $\text{Score}({\bm{y}},{\bm{x}})$ is computed by summing up all the transition scores
and emission scores
where $i$ ranges over all the positions in sentence ${\bm{x}}$, and the substructure space ${\mathbb{U}}({\bm{x}})$ contains all possible position-specific labels $\{y_i\}$ and label pairs $\{(y_{i-1},y_i)\}$.

\subsection{Knowledge Distillation}

\label{sec:kd}
Knowledge distillation is a technique that trains a small student model by encouraging it to imitate the output probability distribution of a large teacher model. The typical KD objective function is the cross-entropy between the output distributions predicted by the teacher model and the student model:
\begin{align}
\mathcal{L}_{\text{KD}}
 = -\sum_{{\bm{y}} \in {\mathbb{Y}}({\bm{x}})} P_t({\bm{y}}|{\bm{x}})\log P_s({\bm{y}}|{\bm{x}}) \label{eq:kd}
\end{align}
where $P_t$ and $P_s$ are the teacher's and the student's distributions respectively.

During training, the student jointly learns from the gold targets and the distributions predicted by the teacher by optimizing the following objective function:
\begin{align}
    \mathcal{L}_{\text{student}} = \lambda \mathcal{L}_\text{KD} + (1-\lambda) \mathcal{L}_\text{target} \nonumber
\end{align}
where $\lambda$ is an interpolation coefficient between the target loss $\mathcal{L}_\text{target}$ and the structural KD loss $\mathcal{L}_\text{KD}$. Following \citet{clark-etal-2019-bam,wang-etal-2020-structure}, one may apply teacher annealing in training by decreasing $\lambda$ linearly from 1 to 0.
Because KD does not require gold labels, unlabeled data can also be used in the KD loss.

\section{Structural Knowledge Distillation}
\label{sec:model}
When performing knowledge distillation on structured prediction, a major challenge is that the structured output space is exponential in size, leading to intractable computation of the KD objective in Eq. \ref{eq:kd}.
However, if the scoring function of the student model can be factorized into scores of substructures (Eq. \ref{eq:2.1}), then we can derive the following factorized form of the structural KD objective.
\begin{align}
&\mathcal{L}_{\text{KD}}
= -\sum_{{\bm{y}} \in {\mathbb{Y}}({\bm{x}})} P_t({\bm{y}}|{\bm{x}}){\log} P_s({\bm{y}}|{\bm{x}}) \nonumber \\
& {=} {-}\smashoperator{\sum_{{\bm{y}} \in {\mathbb{Y}}({\bm{x}})}} P_t({\bm{y}}|{\bm{x}}) \sum_{{\bm{u}} \in {\bm{y}}} \text{Score}_s({\bm{u}},{\bm{x}}) {+} {\log} \mathcal{Z}_s({\bm{x}}) \nonumber\\
&{=} {-} \smashoperator{\sum_{{\bm{y}} \in {\mathbb{Y}}({\bm{x}})}} P_t({\bm{y}}|{\bm{x}}) \smashoperator{\sum_{{\bm{u}} \in {\mathbb{U}}_s({\bm{x}})}}  \bm{1}_{{\bm{u}} \in {\bm{y}}} \text{Score}_s({\bm{u}},{\bm{x}}) {+} {\log} \mathcal{Z}_s({\bm{x}}) \nonumber \\
& {=} {-}\smashoperator{\ssum_{\substack{{\bm{u}} \in {\mathbb{U}}_s({\bm{x}}), {\bm{y}} \in {\mathbb{Y}}({\bm{x}})}}} P_t({\bm{y}}|{\bm{x}})  \bm{1}_{{\bm{u}} \in {\bm{y}}} \text{Score}_s({\bm{u}},{\bm{x}})  {+} {\log} \mathcal{Z}_s({\bm{x}}) \nonumber \\
& {=} {-}\smashoperator{\sum_{{\bm{u}} \in {\mathbb{U}}_s({\bm{x}})}}P_t({\bm{u}}|{\bm{x}}) \text{Score}_s({\bm{u}},{\bm{x}}){+}{\log} \mathcal{Z}_s({\bm{x}})
\label{eq:kl}
\end{align}
where $\bm{1}_\mathrm{condition}$ is 1 if the condition is true and 0 otherwise. From Eq. \ref{eq:kl}, we see that if
${\mathbb{U}}_s({\bm{x}})$ is polynomial in size and $P_t({\bm{u}}|{\bm{x}})$ can be tractably estimated, then the structural KD objective can be tractably computed and optimized.
In the rest of this section, we will show that this is indeed the case for some of the most widely used models in sequence labeling and dependency parsing, two representative structured prediction tasks in NLP. Based on the difference in score factorization between the teacher and student models, we divide our discussion into four scenarios.

\subsection{Teacher and Student Share the Same Factorization Form}

\paragraph{Case 1a: Linear-Chain CRF $\Rightarrow$ Linear-Chain CRF}
In this case, both the teacher and the student are linear-chain CRF models.
An example application is to compress a state-of-the-art CRF model for named entity recognition (NER) that is based on large pretrained contextualized embeddings to a smaller CRF model with static embeddings that is more suitable for fast online serving.

For a CRF student model described in section \ref{sec:str_pred}, if we absorb the emission score $\text{S}_\text{e}(y_i,{\bm{x}})$ into the transition score $\text{S}_\text{t}((y_{i-1},y_i),{\bm{x}})$ at each position $i$, then the substructure space ${\mathbb{U}}_s({\bm{x}})$ contains every two adjacent labels $\{(y_{i-1},y_i)\}$ for $i{=}1,\ldots,n$, with $n$ being the sequence length, and the substructure score is defined as $\text{Score}((y_{i-1},y_i),{\bm{x}})= \text{S}_\text{t}((y_{i-1},y_i),{\bm{x}})+\text{S}_\text{e}(y_i,{\bm{x}})$. The substructure marginal $P_t((y_{i-1}, y_i)|{\bm{x}})$ of the teacher model can be computed by:
\begin{align}

     P_t((y_{i-1}, y_i)|{\bm{x}}) &\propto \alpha(y_{i-1}) \times \beta(y_i)\nonumber\\
     &\times \exp(\text{Score}((y_{i-1},y_i),{\bm{x}}))
     \label{eq:crf2crf}
\end{align}
where $\alpha(y_{i-1})$ and $\beta(y_i)$ are forward and backward scores that can be tractably calculated using the classical forward-backward algorithm.

Comparing with the Posterior KD and Top-K KD of linear-chain CRFs proposed by \citet{wang-etal-2020-structure}, our approach calculates and optimizes the KD objective exactly, while their two KD approaches perform KD either heuristically or approximately. At the formulation level, our approach is based on the marginal distributions of two adjacent labels, while the Posterior KD is based on the marginal distributions of a single label.

\paragraph{Case 1b: Graph-based Dependency Parsing $\Rightarrow$ Dependency Parsing as Sequence Labeling}
In this case, we use the biaffine parser proposed by \citet{dozat-etal-2017-stanfords} as the teacher and the sequence labeling approach proposed by \citet{strzyz-etal-2019-viable} as the student for the dependency parsing task.
The biaffine parser is one of the state-of-the-art models, while the sequence labeling parser provides a good speed-accuracy tradeoff. There is a big gap in accuracy between the two models and therefore KD can be used to improve the accuracy of the sequence labeling parser.

Here we follow the head-selection formulation of dependency parsing without the tree constraint. The dependency parse tree ${\bm{y}}$ is represented by $\langle y_1, \ldots, y_n \rangle$, where $n$ is the sentence length and $y_i=(h_i,l_i)$ denotes the dependency head of the $i$-th token of the input sentence, with $h_i$ being the index of the head token and $l_i$ being the dependency label.
The biaffine parser predicts the dependency head for each token independently. It models separately the probability distribution of the head index $P_t(h_i|{\bm{x}})$ and the probability distribution of the label $P_t( l_i |{\bm{x}})$.
The sequence labeling parser is a MaxEnt model that also predicts the head of each token independently.
It computes $\text{Score}((h_i, l_i), {\bm{x}})$ for each token and applies a softmax function to produce the distribution $P_s((h_i,l_i)|{\bm{x}})$.

Therefore, these two models share the same factorization in which each substructure is a dependency arc specified by $y_i$. ${\mathbb{U}}_s({\bm{x}})$ thus contains all possible dependency arcs among tokens of the input sentence ${\bm{x}}$. The substructure marginal predicted by the teacher can be easily derived as:
\begin{equation}
    P_t((h_i,l_i)|{\bm{x}}) = P_t(h_i|{\bm{x}}) \times P_t(l_i|{\bm{x}})
    \label{eq:dp2sl}
\end{equation}
Note that in this case, the sequence labeling parser uses a MaxEnt decoder, which is locally normalized for each substructure. Therefore, the structural KD objective in Eq. \ref{eq:kl} can be reduced to the following form without the need for calculating the student partition function $\mathcal{Z}_s({\bm{x}})$.
\begin{align}
    \mathcal{L}_{\text{KD}}

    &= -\sum_{{\bm{u}} \in {\mathbb{U}}_s({\bm{x}})}P_t({\bm{u}}|{\bm{x}}) \times \text{log} P_s({\bm{u}}|{\bm{x}})
\end{align}
In all the cases except \textbf{Case 1a} and \textbf{Case 3},
the student model is locally normalized and hence we can follow this form of objective.

\subsection{Student Factorization Produces More Fine-grained Substructures than Teacher Factorization}

\paragraph{Case 2a: Linear-Chain CRF $\Rightarrow$ MaxEnt}
\label{subsubsec:crf2soft}
In this case, we use a linear-chain CRF model as the teacher and a MaxEnt model as the student. Previous work \citep{yang-etal-2018-design,wang-etal-2020-structure} shows that a linear-chain CRF decoder often leads to better performance than a MaxEnt decoder for many sequence labeling tasks. Still, the simplicity and efficiency of the MaxEnt model is desirable. Therefore, it makes sense to perform KD from a linear-chain CRF to a MaxEnt model.

As mentioned in \textbf{Case 1a}, the substructures of a linear-chain CRF model are consecutive labels $\{(y_{i-1},y_i)\}$. In contrast, a MaxEnt model predicts the label probability distribution $P_s(y_i|{\bm{x}})$ of each token independently and
hence the substructure space ${\mathbb{U}}_s({\bm{x}})$ consists of every individual label $\{y_i\}$.
To calculate the substructure marginal of the teacher $P_t(y_i|{\bm{x}})$, we can again utilize the forward-backward algorithm:
\begin{align}
    P_t(y_i|{\bm{x}}) \propto \alpha(y_i) \times \beta(y_i)
\end{align}
where $\alpha(y_{i})$ and $\beta(y_i)$ are forward and backward scores.

\paragraph{Case 2b: Second-Order Dependency Parsing $\Rightarrow$ Dependency Parsing as Sequence Labeling}

The biaffine parser is a first-order dependency parser, which scores each dependency arc in a parse tree independently. A second-order dependency parser scores pairs of dependency arcs with a shared token. The substructures of second-order parsing are therefore all the dependency arc pairs with a shared token.
It has been found that second-order extensions of the biaffine parser often have higher parsing accuracy \citep{wang-etal-2019-second,zhang-etal-2020-efficient,wang-etal-2020-enhanced,wang-tu-2020-second}.
Therefore, we may take a second-order dependency parser as the teacher to improve a sequence labeling parser.

Here we consider the second-order dependency parser of \citet{wang-tu-2020-second}. It employs mean field variational inference to estimate the probabilities of arc existence $P_t(h_i|{\bm{x}})$ and uses a first-order biaffine model to estimate the probabilities of arc labels $P_t(l_i|{\bm{x}})$.
Therefore, the substructure marginal can be calculated in the same way as Eq. \ref{eq:dp2sl}.

\subsection{Teacher Factorization Produces More Fine-grained Substructures than Student Factorization}
\label{subsec:soft2crf}
\paragraph{Case 3: MaxEnt $\Rightarrow$ Linear-Chain CRF}
Here we consider KD in the opposite direction of \textbf{Case 2a}. An example application is zero-shot cross-lingual NER.
Previous work \citep{pires-etal-2019-multilingual,wu-dredze-2019-beto} has shown that multilingual BERT (M-BERT) has strong zero-shot cross-lingual transferability in NER tasks. Many such models employ a MaxEnt decoder.
In scenarios requiring fast speed and low computation cost, however, we may want to distill knowledge from such models to a model with much cheaper static monolingual embeddings while compensating the performance loss with a linear-chain CRF decoder.

As described in \textbf{Case 1a}, the substructures of a linear-chain CRF model are consecutive labels $\{(y_{i-1},y_i)\}$. Because of the label independence and local normalization in the MaxEnt model, the substructure marginal of the MaxEnt teacher is calculated by:
\begin{align}
    P_t((y_{i-1}, y_i)|{\bm{x}}) = P_t(y_{i-1}|{\bm{x}})P_t(y_i|{\bm{x}})
    \label{max2crf}
\end{align}

\subsection{Factorization Forms From Teacher and Student are Incompatible}
\paragraph{Case 4: NER as Parsing $\Rightarrow$ MaxEnt}
Very recently, \citet{yu-etal-2020-named} propose to solve the NER task as graph-based dependency parsing and achieve state-of-the-art performance. They represent each named entity with a dependency arc from the first token to the last token of the named entity, and represent the entity type with the arc label. However, for the flat NER task (i.e., there is no overlapping between entity spans), the time complexity of this method is higher than commonly used sequence labeling NER methods. In this case, we take a parsing-based NER model as our teacher and a MaxEnt model with the BIOES label scheme as our student.

The two models adopt very different representations of NER output structures. The parsing-based teacher model represents an NER output of a sentence with a set of labeled dependency arcs and defines its score as the sum of arc scores. The MaxEnt model represents an NER output of a sentence with a sequence of BIOES labels and defines its score as the sum of token-wise label scores. Therefore, the factorization forms of these two models are incompatible.

Computing the substructure marginal of the teacher $P_t(y_i|{\bm{x}})$, where $y_i \in \{B_l, I_l, E_l, S_l, O | l \in L\}$ and $L$ is the set of entity types, is much more complicated than in the previous cases. Take $y_i = B_l$ for example. $P_t(y_i = B_l |{\bm{x}})$ represents the probability of the $i$-th word being the beginning of a multi-word entity of type `$l$’. In the parsing-based teacher model, this probability is proportional to the summation of exponentiated scores of all the output structures that contain a dependency arc of label `$l$’ with the $i$-th word as its head and with its length larger than 1. It is intractable to compute such marginal probabilities by enumerating all the output structures, but
we can tractably compute them using dynamic programming. See supplementary material for a detailed description of our dynamic programming method.

\section{Experiments}
\label{sec:exp}
We evaluate our approaches described in Section \ref{sec:model} on NER (\textbf {Case 1a, 2a, 3, 4}) and dependency parsing (\textbf{Case 1b, 2b}).

\subsection{Settings}
\label{sec:set}
\paragraph{Datasets} We use CoNLL 2002/2003 datasets \citep{tjong-kim-sang-2002-introduction,tjong-kim-sang-de-meulder-2003-introduction} for {\bf Case 1a}, \textbf{2a} and \textbf{4}, and use WikiAnn datasets \citep{pan-etal-2017-cross} for  {\bf Case 1a}, \textbf{2a}, \textbf{3}, and \textbf{4}. The CoNLL datasets contain the corpora of four Indo-European languages. We use the same four languages from the WikiAnn datasets. For cross-lingual transfer in \textbf{Case 3}, we use the four Indo-European languages as the source for the teacher model and additionally select four languages from different language families as the target for the student models.\footnote{The four languages from the CoNLL datasets are Dutch, English, German and Spanish and the four target languages for \textbf{Case 3} are Basque, Hebrew, Persian and Tamil. We use ISO 639-1 language codes (\url{https://en.wikipedia.org/wiki/List_of_ISO_639-1_codes}) to represent each language.}

We use the standard training/development/test split for the CoNLL datasets. For WikiAnn, we follow the sampling of \citet{wang-etal-2020-structure} with 12000 sentences for English and 5000 sentences for each of the other languages. We split the datasets by 3:1:1 for training/development/test. For \textbf{Case 1b} and \textbf{2b}, we use Penn Treebank (PTB) 3.0 and follow the same pre-processing pipeline as in \citet{ma-etal-2018-stack}.
For unlabeled data, we sample sentences that belong to the same languages of the labeled data from the WikiAnn datasets for {\bf Case 1a}, \textbf{2a} and {\bf 4} and we sample sentences from the target languages of WikiAnn datasets for \textbf{Case 3}. We use the BLLIP corpus\footnote{Brown Laboratory for Linguistic Information Processing (BLLIP) 1987-89 WSJ Corpus Release 1. } as the unlabeled data for \textbf{Case 1b} and \textbf{2b}.

\paragraph{Models}

For the student models in all the cases, we use fastText \citep{bojanowski2017enriching} word embeddings and character embeddings as the word representation. For {\bf Case 1a}, \textbf{2a} and {\bf 4}, we concatenate the multilingual BERT, Flair \citep{akbik-etal-2018-contextual}, fastText embeddings and character embeddings \citep{santos2014learning} as the word representations for stronger monolingual teacher models \cite{wang-etal-2020-more}. For {\bf Case 3}, we use M-BERT embeddings for the teacher. Also for \textbf{Case 3}, we fine-tune the teacher model on the training set of the four Indo-European languages from the WikiAnn dataset and train student models on the four additional languages. For the teacher models in {\bf Case 1b} and {\bf 2b}, we simply use the same embeddings as the student because there is already huge performance gap between the teacher and student in these settings and hence we do not need strong embeddings for the teacher to demonstrate the utility of KD.

\paragraph{Baselines}
We compare our Structural KD (\textit{Struct. KD}) with training without KD (\textit{w/o KD}) as well as existing KD approaches.
In \textbf{Case 1a}, the \textit{Pos. KD} baseline is the Posterior KD approach for linear-chain CRFs proposed by \citet{wang-etal-2020-structure}. They also propose \textit{Top-K KD} but have shown that it is inferior to \textit{Pos. KD}. For experiments using unlabeled data in all the cases, in addition to labeled data, we use the teacher's prediction on the unlabeled data as pseudo labeled data to train the student models. This can be seen as the \textit{Top-1 KD} method\footnote{We do not predict pseudo labels for the labeled data, because we find that the teacher models' predictions on the labeled training data have approximately 100\% accuracy in most of the cases. }.
In Case \textbf{2a} and \textbf{3}, where we perform KD between CRF and MaxEnt models, we run a reference baseline that replaces the CRF teacher or student model with a MaxEnt model and performs token-level KD (\textit{Token KD}) of MaxEnt models that optimizes the cross entropy between the teacher and student label distributions at each position.

\begin{table}[t]
\small
\centering
\setlength\tabcolsep{1pt}
\begin{tabular}{l|cc|c|cc|c|cc}
\hlineB{4}

Case     & \multicolumn{2}{c|}{\textbf{1a}} & \textbf{1b} & \multicolumn{2}{c|}{\textbf{2a}} & \textbf{2b} & \multicolumn{2}{c}{\textbf{4}} \\

Labeled      & CoN            & Wiki            & PTB              & CoN            & Wiki            & PTB              & CoN           & Wiki           \\ \hline\hline
Teacher      & 89.15            & 88.52              & 95.96            & 89.15            & 88.52              & 96.04            & 88.57           & 88.38             \\
\hline
w/o KD        & 84.70            & 83.31              & 89.85            & 83.87            & 80.86              & 89.85            & 83.87           & 80.86             \\
Pos. KD & 85.27            & 83.73              & -                & -                & -                  & -                & -               & -                 \\
Struct. KD    & \textbf{85.35}   & \textbf{84.12}     & \textbf{91.83}   & \textbf{84.50}   & \textbf{82.23}     & \textbf{91.78}   & \textbf{84.28}  & \textbf{81.45}    \\
\hlineB{4}
\end{tabular}
\caption{Averaged F1 scores for NER and labeled attachment scores (LAS) for dependency parsing on labeled datasets. CoN: CoNLL datasets.}
\label{tab:labeled}
\end{table}

\begin{table}[t]
\small
\centering
\setlength\tabcolsep{1.5pt}
\begin{tabular}{l|c|c|c|c|c|c}
\hlineB{4}
Case       & \textbf{1a} & \textbf{1b} & \textbf{2a} & \textbf{2b} & \textbf{3} & \textbf{4} \\
Labeled+Unlabeled        & Wiki           & PTB           & Wiki           & PTB            & Wiki U          & Wiki          \\
\hline\hline
Teacher        & 88.52            & 95.96            & 88.52            & 96.04            & 56.01           & 88.38           \\
\hline
Top-1       & 84.19            & 90.03            & 82.40             & 90.03            & 41.11           & 82.10            \\
Pos. KD + Top-1   & 84.91            & -                & -                & -                & -               & -               \\
Struct. KD + Top-1 & \textbf{85.24}   & \textbf{91.98}   & \textbf{85.24}   & \textbf{91.94}   & \textbf{45.28}  & \textbf{82.44}  \\ \hlineB{4}
\end{tabular}
\caption{Average F1 score of NER and labeled attachment scores (LAS) for dependency parsing with both labeled and unlabeled data. Wiki U means that the training data of this case contains only the unlabeled data. }
\label{tab:unlabel}
\end{table}

\paragraph{Training}
For MaxEnt and linear-chain CRF models, we use the same hyper-parameters as in \citet{akbik-etal-2018-contextual}. For dependency parsing, we use the same hyper-parameters as in \citet{wang-tu-2020-second} for teacher models and \citet{strzyz-etal-2019-viable} for student models. For M-BERT fine-tuning in {\bf Case 3}, we mix the training data of the four source datasets and train the teacher model with the AdamW optimizer \citep{loshchilov2018decoupled} with a learning rate of $5{\times} 10^{-5}$ for 10 epochs. We tune the KD temperature in $\{1,2,3,4,5\}$ and the loss interpolation annealing rate in $\{0.5,1.0,1.5\}$.
For all experiments, we train the models for 5 runs with a fixed random seed for each run.

\subsection{Results}
\label{sec:results}

\begin{table}[t]
\small
\centering
\setlength\tabcolsep{2pt}
\begin{tabular}{l||ccc|c}
\hlineB{4}
                 & \multicolumn{3}{c|}{\textbf{Case 2a}}                    & \textbf{Case 3}         \\
                 & CoNLL         & WikiAnn        & Wiki+U         & Wiki U         \\ \hline\hline
MaxEnt Teacher & 88.65         & 87.41          & 87.41          & 56.01          \\
CRF Teacher      & 89.15         & 88.52          & 88.52          & -            \\
\hline
Token. KD        & 84.25         & 82.09          & 83.07          & 38.42          \\
Struct. KD       & \textbf{84.50} & \textbf{82.23} & \textbf{83.34} & \textbf{45.28} \\ \hlineB{4}
\end{tabular}
\caption{Comparing with reference baselines on NER task. Wiki+U means the training data comprises labeled and unlabeled WikiAnn data and Wiki U means that the training data of this case contains only the unlabeled data.}
\label{tab:baseline_ref}
\end{table}

Table \ref{tab:labeled} shows the experimental results with labeled data only and \ref{tab:unlabel} shows the experimental results with additional 3000 unlabeled sentences.
The results show that our structural KD approaches outperform the baselines in all the cases. Table \ref{tab:baseline_ref} compares \textit{Struct. KD} with \textit{Token KD}, the reference baseline based on MaxEnt models.
For \textbf{Case 2a}, which involves a MaxEnt student, \textit{Struct. KD} with a CRF teacher achieves better results than \textit{Token KD} with a MaxEnt teacher. For \textbf{Case 3}, which involves a MaxEnt teacher, \textit{Struct. KD} with a CRF student achieves better results than \textit{Token KD} with a MaxEnt student.
These results are to be expected because \textit{Struct. KD} makes it possible to apply exact knowledge distillation with a more capable teacher or student.
In all the experiments, we run Almost Stochastic Dominance proposed by \citet{dror-etal-2019-deep} with a significance level of $0.05$ and find that the advantages of our structural KD approaches are significant. Please refer to Appendix for more detailed results.

\begin{table}[t]
\setlength\tabcolsep{5pt}
\centering
\small
\begin{tabular}{l||cccc|c}
\hlineB{4}
&	 de    	 &   en &	  es    	 & nl    &	  Avg.  \\
\hline
\multicolumn{6}{c}{\textbf{Case 1a}}\\
\hline
w/o KD$^\dagger$         &      	82.16&	90.13&	88.06&	89.11&	 87.36\\
Top-WK KD$^\dagger$  &    	82.15&	90.52&	88.64&	89.24&	   87.64\\
Pos. KD$^\dagger$  &            	82.22&	90.68&	88.57&	89.41&	 87.72\\
Pos.+Top-WK$^\dagger$ &	 \textbf{82.31} &	90.53	&88.66&	89.58&	 87.77\\
Struct. KD& 	82.28&	 \textbf{90.86}& 	 \textbf{88.67}& 	 \textbf{90.07}& 	 \textbf{87.97} \\
\hline
\multicolumn{6}{c}{\textbf{Case 2a}}\\
\hline
w/o KD$^\dagger$ &	 \textbf{81.40}  &	90.08 &	87.72 &	88.99	 & 87.05 \\
Token KD$^\dagger$  &	81.30 &	90.02 &	88.24 &	88.87 &	 87.11 \\
Struct. KD  &	81.27 &	 \textbf{90.25}  &	 \textbf{88.64}  &	 \textbf{89.14}  &	 \textbf{87.32}\\
\hlineB{4}
\end{tabular}
\caption{A comparison of KD approaches for multilingual NER. $\dagger$: Results are from \citet{wang-etal-2020-structure}.}
\label{tab:ner}
\end{table}

\subsection{Multilingual NER Experiments}
There is a recent increase of interest in training multilingual NER models \citep{tsai-etal-2019-small,mukherjee-hassan-awadallah-2020-xtremedistil} because of the strong generalizability of M-BERT on multiple languages. Existing work explored knowledge distillation approaches to train fast and effective multilingual NER models with the help of monolingual teachers \citep{wang-etal-2020-structure}. To show the effectiveness of structural KD in the multilingual NER setting, we compare our approaches with those reported by \citet{wang-etal-2020-structure}.
Specifically, the monolingual teachers are always CRF models, and the multilingual student is either a CRF model (\textbf{Case 1a}) or a MaxEnt model (\textbf{Case 2a}). \citet{wang-etal-2020-structure} report results of the \textit{Top-WK KD} (a weighted version of \textit{Top-K KD}) and \textit{Pos. KD} approaches for \textbf{Case 1a} and the reference baseline \textit{Token KD} (with a MaxEnt teacher) for \textbf{Case 2a}. We follow their experimental settings when running our approach.

The experimental results in Table \ref{tab:ner} show the effectiveness of \textit{Struct. KD} in both cases.
In \textbf{Case 1a}, our approach is stronger than both \textit{Top-WK KD} and \textit{Pos. KD} as well as the mixture of the two approaches on average.
In \textbf{Case 2a}, \textit{Struct. KD}
not only outperforms \textit{Token KD}, but also makes the MaxEnt student competitive with the CRF student without KD (87.32 vs. 87.36).

\begin{filecontents}{zs_maxent2crf.dat}
unlabel	MaxEnt	CRF	MaxEnt->MaxEnt	Structural Teacher baseline
3	37.0515	41.1075	38.417	45.2785 56.01 37.0515
10	44.0215	48.457	45.8685	51.678 56.01 37.0515
30	53.794	57.2435	54.827	59.2375 56.01 37.0515
\end{filecontents}

\begin{filecontents}{graph2seq_parsing.dat}
unlabel	Baseline	Structural Teacher baseline
0	89.854	91.828 94.24 89.854
3	90.028	91.976 94.24 89.854
10	90.622	92.236 94.24 89.854
30	91.44	92.632 94.24 89.854
50	91.9925	92.996 94.24 89.854
100	92.61	93.282 94.24 89.854
\end{filecontents}

\begin{filecontents}{crf2crf.dat}
unlabel	crf	posterior	structural teacher baseline
0	83.31	83.7275	84.1225 88.5225 83.31
3	84.185	85.005	85.2075 88.5225 83.31
10	84.9	85.4525	85.6625 88.5225 83.31
30	86.24	86.5075	86.64 88.5225 83.31
\end{filecontents}

\begin{figure*}[!ht]
\centering
\begin{tikzpicture}

    \node at (-4.7,-0.65) {(a) Case 1a};
    \node at (0.7,-0.65) {(b) Case 1b};
    \node at (5.9,-0.65) {(c) Case 3};

    \begin{axis}[
        xshift=-7cm,

        name=ner,

        width=0.38\textwidth,
        height=0.25\textwidth,

        legend columns=2,
        legend pos=north west,
        legend style={font=\tiny,yshift=-0.3cm},
        tick label style={font=\small},

        xtick={0,3,10,30},
        ylabel style={yshift=-0.5cm},
        ]
        \addplot[red,mark=square] table[x=unlabel,y=crf] {crf2crf.dat};
        \addplot[black,mark=*] table[x=unlabel,y=posterior] {crf2crf.dat};
        \addplot[blue,mark=square*] table[x=unlabel,y=structural] {crf2crf.dat};
        \addplot[black,dashed,line width=0.4mm] table[x=unlabel,y=teacher] {crf2crf.dat};
        \addplot[black,dotted,line width=0.4mm] table[x=unlabel,y=baseline] {crf2crf.dat};
        \legend{Top-1 KD, Pos. KD, Struct. KD}
    \end{axis}
    \begin{axis}[
        xshift=5.3cm,

        at={(ner.south west)},
        name=zeroshot,

        width=0.38\textwidth,
        height=0.25\textwidth,

        legend columns=1,
        legend pos=south east,
        legend style={font=\tiny,yshift=0.3cm},
        tick label style={font=\small},

         xtick={0,10,30,50,100},

        ylabel style={yshift=-0.5cm},
        ]
        \addplot[red,mark=square] table[x=unlabel,y=Baseline] {graph2seq_parsing.dat};
        \addplot[blue,mark=square*] table[x=unlabel,y=Structural] {graph2seq_parsing.dat};
        \addplot[black,dashed,line width=0.4mm] table[x=unlabel,y=Teacher] {graph2seq_parsing.dat};
        \addplot[black,dotted,line width=0.4mm] table[x=unlabel,y=baseline] {graph2seq_parsing.dat};
        \legend{Top-1 KD, Struct. KD}
    \end{axis}
    \begin{axis}[

        at={(zeroshot.south west)},
        xshift=5.3cm,

        width=0.38\textwidth,
        height=0.25\textwidth,

        legend columns=1,
        legend pos=south east,
        legend style={font=\tiny},
        tick label style={font=\small},

        xtick={0, 3,10,30},
        ylabel style={yshift=-0.5cm},
        ]

        \addplot[red,mark=square] table[x=unlabel,y=CRF] {zs_maxent2crf.dat};
        \addplot[black,mark=*] table[x=unlabel,y=MaxEnt->MaxEnt] {zs_maxent2crf.dat};
        \addplot[blue,mark=square*] table[x=unlabel,y=Structural] {zs_maxent2crf.dat};
        \addplot[black,dashed,line width=0.4mm] table[x=unlabel,y=Teacher] {zs_maxent2crf.dat};

        \legend{Top-1 KD, Token KD, Struct. KD}

    \end{axis}
\end{tikzpicture}
\caption{The accuracy of structural KD and the baselines on different amounts of unlabeled data in three cases. The x-axis represents the amount of unlabeled data in thousand and the y-axis represents the accuracy. The dashed lines are the accuracy of the teacher models. The dotted lines are the accuracy of the baseline models without any knowledge from the teachers.}
\label{fig:performance_curve}
\end{figure*}

\section{Analysis}
\subsection{Amount of Unlabeled Data}

We compare our approaches with the baselines with different amounts of unlabeled data for \textbf{Case 1a}, \textbf{1b} and \textbf{3}, which are cases that apply in-domain unlabeled data for NER and dependency parsing, and cross-lingual unlabeled data for NER. We experiment with more unlabeled data for \textbf{Case 1b} than for the other two cases because the labeled training data of PTB is more than $10$ times larger than the labeled NER training data in \textbf{Case 1a} and \textbf{3}. Results are shown in Figure \ref{fig:performance_curve}.
The experimental results show that our approaches consistently outperform the baselines, though the performance gaps between them become smaller when the amount of unlabeled data increases. Comparing the performance of the students with the teachers, we can see that in \textbf{Case 1a} and \textbf{1b}, the gap between the teacher and the student remains large even with the largest amount of unlabeled data. This is unsurprising considering the difference in model capacity between the teacher and the student. In \textbf{Case 3}, however, we find that when using 30,000 unlabeled sentences, the CRF student models can even outperform the MaxEnt teacher model, which shows the effectiveness of CRF models on NER.

\subsection{Temperature in Structural Knowledge Distillation}
A frequently used KD technique is dividing the logits of probability distributions of both the teacher and the student by a temperature in the KD objective \citep{44873}. Using a higher temperature produces softer probability distributions and often results in higher KD accuracy. In structural KD, there are two approaches to applying the temperature to the teacher model, either globally to the logit of $P_t({\bm{y}}|{\bm{x}})$ (i.e., $\text{Score}_t({\bm{y}},{\bm{x}})$) of the full structure ${\bm{y}}$, or locally to the logit of $P_t({\bm{u}}|{\bm{x}})$ of each student substructure ${\bm{u}}$.
We empirically compare these two approaches in \textbf{Case 1a} with the same setting as in Section \ref{sec:set}. Table \ref{tab:temperature} shows that the local approach results in better accuracy for all the languages. Therefore, we use the local approach by default in all the experiments.

\begin{table}[!t]

\small
\centering
\begin{tabular}{l||cccc|c}
\hlineB{4}
 & de    & en    & es    & nl    & Avg.   \\
\hline\hline
CRF    & 75.37 & 91.21 & 86.55 & 85.67 & 84.70 \\
Global      & 75.67 & 91.11 & 86.72 & 85.92 & 84.85 \\
Local       & \textbf{76.61} & \textbf{91.41} & \textbf{87.20} & \textbf{86.19} & \textbf{85.35}\\
\hlineB{4}
\end{tabular}
\captionof{table}{Comparison of the global and local temperature application approaches on CoNLL NER.}
\label{tab:temperature}
\end{table}
\begin{table}[!t]
\small \centering
\begin{tabular}{l||cc}
\hlineB{4}
                  & CoNLL   & WikiAnn \\
\hline\hline
CRF                & \textbf{89.15}  & \textbf{88.52} \\
CRF-Mrg.      & 89.08 & 88.41 \\
NER-Par.           & 88.57   & 88.38 \\
NER-Par.-Mrg. & 87.40   & 86.82 \\
MaxEnt            & 88.65 & 87.41  \\
\hlineB{4}
\end{tabular}
\captionof{table}{Averaged F1 score of teachers and it's marginal distributions. -Mrg.: Marginal distribution, NER-Par.: NER as parsing \citep{yu-etal-2020-named}.}
\label{tab:teachers}
\end{table}

\subsection{Comparison of Teachers}
In \textbf{Case 2a} and \textbf{Case 4}, we use the same MaxEnt student model but different types of teacher models. Our structural KD approaches in both cases compute the marginal distribution $P_t(y_i|{\bm{x}})$ of the teacher at each position $i$ following the substructures of the MaxEnt student, which is then used to train the student substructure scores. We can evaluate the quality of the marginal distributions by taking their modes as label predictions and evaluating their accuracy.
In Table \ref{tab:teachers}, we compare the accuracy of the CRF teacher and its marginal distributions from \textbf{Case 2a}, the NER-as-parsing teacher and its marginal distributions from \textbf{Case 4}, and the MaxEnt teacher which is the KD baseline in \textbf{Case 2a}.
First, we observe that for both CRF and NER-as-parsing, predicting labels from the marginal distributions leads to lower accuracy. This is to be expected because such predictions do not take into account correlations between adjacent labels. While predictions from marginal distributions of the CRF teacher still outperform MaxEnt, those of the NER-as-parsing teacher clearly underperform MaxEnt. This provides an explanation as to why \textit{Struct. KD} in \textbf{Case 4} has equal or even lower accuracy than the \textit{Token KD} baseline in \textbf{Case 2a} in Table \ref{tab:baseline_ref}.

\section{Related Work}

\subsection{Structured Prediction}

In this paper, we use sequence labeling and dependency parsing as two example structured prediction tasks.
In sequence labeling, a lot of work applied the linear-chain CRF and achieved state-of-the-art performance in various tasks \citep{ma-hovy-2016-end,akbik-etal-2018-contextual,liu-etal-2019-gcdt,yu-etal-2020-named,wei-etal-2020-dont,wang2020automated,wang2021improving}. Meanwhile, a lot of other work used the MaxEnt layer instead of the CRF for sequence labeling \citep{devlin-etal-2019-bert,conneau-etal-2020-unsupervised,wang-etal-2020-ain} because MaxEnt makes it easier to fine-tune pretrained contextual embeddings in training. Another advantage of MaxEnt in comparison with CRF is its speed. \citet{yang-etal-2018-design} showed that models equipped with the CRF are about two times slower than models with the MaxEnt layer in sequence labeling.
In dependency parsing, recent work shows that second-order CRF parsers achieve significantly higher accuracy than first-order parsers \citep{wang-etal-2019-second,zhang-etal-2020-efficient}. However, the inference speed of second-order parsers is much slower. \citet{zhang-etal-2020-efficient} showed that second-order parsing is four times slower than the simple head-selection first-order approach \citep{dozat2016deep}.
Such speed-accuracy tradeoff as seen in sequence labeling and dependency parsing also occurs in many other structured prediction tasks. This makes KD an interesting and very useful technique that can be used to circumvent this tradeoff to some extent.

\subsection{Knowledge Distillation in Structured Prediction}
KD has been applied in many structured prediction tasks in the fields of NLP, speech recognition and computer vision, with applications such as neural machine translation \citep{kim-rush-2016-sequence,tan2018multilingual},  sequence labeling \citep{tu-gimpel-2019-benchmarking,wang-etal-2020-structure},
connectionist temporal classification \citep{Huang2018},
image semantic segmentation \citep{liu2019structured} and so on.
In KD for structured prediction tasks, how to handle the exponential number of structured outputs is a main challenge. To address this difficult problem, recent work resorts to approximation of the KD objective. \citet{kim-rush-2016-sequence} proposed sequence-level distillation through predicting K-best sequences of the teacher in neural machine translation. \citet{kuncoro-etal-2016-distilling} proposed to use multiple greedy parsers as teachers and generate the probability distribution at each position through voting. Very recently, \citet{wang-etal-2020-structure} proposed structure-level knowledge distillation for linear-chain CRF models in multilingual sequence labeling. During the distillation process, teacher models predict the Top-K label sequences as the global structure information or the posterior label distribution at each position as the local structural information, which is then used to train the student. Besides approximate approaches, an alternative way is using models that make local decisions and performing KD on these local decisions. \citet{anderson-gomez-rodriguez-2020-distilling} formulated dependency parsing as a head-selection problem and distilled the distribution of the head node at each position. \citet{tsai-etal-2019-small} proposed MiniBERT through distilling the output distributions of M-BERT models of the MaxEnt classifier. Besides the output distribution, \citet{mukherjee-hassan-awadallah-2020-xtremedistil} further distilled the hidden representations of teachers.

\section{Conclusion}
In this paper, we propose structural knowledge distillation, which transfers knowledge between structured prediction models. We derive a factorized form of the structural KD objective and make it tractable to compute and optimize for many typical choices of teacher and student models. We apply our approach to four KD scenarios with six cases for sequence labeling and dependency parsing. Empirical results show that our approach outperforms baselines without KD as well as previous KD approaches. With sufficient unlabeled data, our approach can even boost the students to outperform the teachers in zero-shot cross-lingual transfer.

\subsubsection*{Acknowledgments}
This work was supported by the National Natural Science Foundation of China (61976139) and by Alibaba Group through Alibaba Innovative Research Program.