# Single band VLBI absolute astrometry

Leonid Petrov NASA Goddard Space Flight Center
Code 61A, 8800 Greenbelt Rd, Greenbelt, 20771 MD, USA Leonid Petrov
<EMAIL_ADDRESS>
(Received 2022, November 09)

###### Abstract

The ionospheric path delay impacts single-band very long baseline
interferometry (VLBI) group delays, which limits their applicability for
absolute astrometry. I consider two important cases: when observations are
made simultaneously at two bands, but delays at only one band are available
for a subset of observations and when observations are made at one band only
by design. I developed optimal procedures of data analysis for both cases
using Global Navigation Satellite System (GNSS) ionosphere maps, provided a
stochastic model that describes ionospheric errors, and evaluated their impact
on source position estimates. I demonstrate that the stochastic model is
accurate at a level of 15%. I found that using GNSS ionospheric maps as is
introduces serious biases in estimates of declinations and I developed the
procedure that almost eliminates them. I found serendipitously that GNSS
ionospheric maps have multiplicative errors and have to be scaled by 0.85 in
order to mitigate the declination bias. A similar scale factor was found in
comparison of the vertical total electron contents from satellite altimetry
against GNSS ionospheric maps. I favor interpretation of this scaling factor
as a manifestation of the inadequacy of the thin shell model. I showed in this
study that we are able to model the ionospheric path delay to the extent that
no systematic errors emerge and we are able to adequately assess the
contribution of the ionosphere-driven random errors on source positions. This
makes single-band absolute astrometry a viable option that can be used for
source position determination.

Astrometry, VLBI, Earth ionosphere

††journal: Astronomical Journal††software: PIMA

##  1 Introduction

Nowadays it became customary when speaking about VLBI astrometry to implicitly
assume differential astrometry (see, for instance, Reid & Honma, 2014; Rioja &
Dodson, 2020). Differential VLBI astrometry is based on analysis of fringe
phase differences between a target and a calibrator located within several
degrees of each other. This approach allows to estimate a position offset
between a target and a calibrator with a precision of 30–50 $\mu$as. This high
precision becomes feasible because the contribution of path delay in the
atmosphere is diluted by the factor equal to a target–calibrator separation
expressed in radians. At frequencies below approximately 5 GHz the
contribution of the ionosphere usually dominates, at higher frequencies its
contribution is less than the contribution of the neutral atmosphere. Since
the number of potential calibrators grows with a decrease of frequency,
observations at low frequencies 1–2 GHz are usually made with shorter
target–calibrator arcs. This better compensates for the residual ionospheric
contribution which is larger at low frequencies. Therefore, accuracy of
differential astrometry is approximately the same within a range of
frequencies 1–43 GHz.

I should stress the fact that is often omitted in literature: differential
astrometry alone cannot provide position of a target. For some applications,
such as parallax determination, measurement of proper motion, or estimation of
an orbital motion of a galactic object, knowledge of a precise position is not
required. But there are applications that do require knowledge of precise
positions. Source position can be determined by observing either a target or a
calibrator with the absolute astrometry method. In a latter case, position of
the target can be derived from the position of the calibrator and the position
target/calibrator offset from differential VLBI observations under assumption
that a) positions of the calibrator are stable between epochs of observations;
b) a position offset of the calibrator at the frequency of differential
observations with respect to its position at the frequency of absolute
astrometry observations is well known.

Method of VLBI absolute astrometry involves observations of many active
galactic nuclei approximately uniformly distributed over the sky. Data
analysis of group delays derived from these observations is usually performed
in the accumulative mode, i.e. all VLBI absolute astrometry and geodesy
observations are processed in a single least square solution for estimation of
source coordinates, orientation parameters, stations positions, and nuisance
parameters, such as atmospheric path delay in zenith direction and clock
function. The errors in source position adjustments are due to the thermal
noise and the inaccuracy of modeling path delay in the atmosphere. They vary
in a wide range from 0.05 to 50 mas depending on source flux density, network
geometry, and the number of observables, with 1 mas error being typical.

Since the contribution of the ionosphere to group delay is reciprocal to the
square of the observing frequency, usually absolute astrometry experiments are
performed at two widely separated bands, 2.3/8.4 or 4.3/7.6 GHz (Petrov,
2021). Processing of dual-band data allows to eliminate effectively the
contribution of the ionosphere. However, there are situations when absolute
astrometry observations are available only at one band.

We know that compared with dual-band observations, single band absolute
astrometry observations are affected by the contribution of the path delay in
the ionosphere. Two questions arise: 1) which data analysis strategy does
provide source coordinate estimates with the lowest uncertainties and 2) how
to account for the contribution of errors in the ionosphere modeling to
reported source coordinate uncertainties? The goal of this study is to provide
answers to these questions. I took a trial dataset of over 4 millions of dual-
band twenty-four hour VLBI observations publicly available at the US National
Radio Astronomy Observatory (NRAO) archive and used it as a testbed. I dropped
existing observations at the second band during data analysis and compared
results of single-band data against the reference solution that used both
bands.

##  2 VLBI dataset used for trials

I used a dataset of 263 twenty four hour experiments at Very Large Baseline
Array (VLBA) network since April 1998 though May 2021. All these data are
publicly available at NRAO archive111https://data.nrao.edu. The dataset
consists of observing sessions under regular VLBI geodesy program RDV (Petrov
et al., 2009), astrometric VCS-II program (Gordon et al., 2016), its follow-
ups VCS-III and VCS-IV, and geodetic CONT17 campaign (Behrend et al., 2020).
The motivation of this choice is to have a long history of observations, a
homogeneous network, and both short and long baselines. VLBI absolute
astrometry programs at declinations above $-40^{\circ}$ are almost exclusively
run with VLBA. Therefore, conclusions made from processing trial runs at the
VLBA can be propagated directly to the past and future astrometry programs at
that network.

##  3 Modeling the ionospheric contribution to path delay

###  3.1 Dual-band observations

The impact of the ionosphere dispersiveness on fringe phase is reciprocal to
frequency in the first approximation. Therefore, fringe phase in channel $i$
in the presence of the ionosphere becomes

$\displaystyle\phi_{i}=2\pi\uptau_{p}\,f_{0}+\uptau_{g}\,(f_{i}-f_{0})-\frac{\displaystyle\alpha}{\displaystyle
f_{i}},$ (1)

where $\uptau_{p}$ and $\uptau_{g}$ are phase and group delays, $f_{i}$ is the
frequency of the $i$th spectral channel, $f_{0}$ is the reference frequency
and

$\displaystyle\alpha=\frac{\displaystyle e^{2}}{\displaystyle
8\,\pi^{2}\,c\,m_{e}\,\epsilon_{o}}\left(\int N_{v}\,ds_{1}-\int
N_{v}\,ds_{2}\right),$ (2)

where $N_{v}$ — electron density, $e$ — charge of an electron, $m_{e}$ — mass
of an electron, $\epsilon_{o}$ — permittivity of free space, and $c$ —
velocity of light in vacuum. Integration is carried along the line of sight.
Having substituted values of constants (Klobuchar, 1996) and expressing the
total electron contents along the line of sight $\int N_{v}\,ds$ in
$\mbox{1}\cdot 10^{16}$ electrons/$m^{2}$ (so-called TEC units or TECU), we
arrive to $\alpha=\mbox{1.345}\cdot 10^{9}\>\>\mbox{sec}$/TECU.

Phase and group delay are computed using fringe phases $\phi_{i}$ with weights
$w_{i}$ using least squares. The result can be expressed analytically after
some algebra:

$\displaystyle\uptau_{\rm gi}=\uptau_{\rm
if}+\frac{\displaystyle\alpha}{\displaystyle f_{e}^{2}}\,\rm TEC,$ (3)

where $\uptau_{\rm if}$ is the ionosphere-free group delay, $\rm TEC$ is $\int
N_{v}\,ds$ expressed in TEC units, and $f_{e}$ is the effective ionospheric
frequency

$\displaystyle
f_{e}=\sqrt{\frac{\displaystyle\displaystyle\sum_{i}^{n}w_{i}\cdot\displaystyle\sum_{i}^{n}w_{i}(f_{i}-f_{0})^{2}\>-\>\left(\displaystyle\sum_{i}^{n}w_{i}(f_{i}-f_{0})\right)^{2}\,}{\displaystyle\displaystyle\sum_{i}^{n}w_{i}(f_{i}-f_{0})\displaystyle\sum_{i}^{n}\frac{w_{i}}{f_{i}}\>-\>\displaystyle\sum_{i}^{n}w_{i}\cdot\displaystyle\sum_{i}^{n}w_{i}\frac{(f_{i}-f_{0})}{f_{i}}}}$
(4)

that depends on weights of spectral channels $w_{i}$. Typically, the effective
ionospheric frequency is within several percents of the central frequency of
the observing band.

The best way to mitigate the impact of the ionosphere on group delay is to
observe simultaneously at two or more widely separated frequency bands. Then
the following linear combination of two group delays at the upper and lower
bands, $\uptau_{u}$ and $\uptau_{l}$, respectively is ionosphere free:

$\displaystyle\uptau_{if}=\frac{\displaystyle f_{u}^{2}}{\displaystyle
f_{u}^{2}-f_{l}^{2}}\uptau_{u}-\frac{\displaystyle f_{u}^{2}}{\displaystyle
f_{u}^{2}-f_{l}^{2}}\uptau_{l}.$ (5)

Here $f_{u}$ and $f_{l}$ are effective ionospheric frequencies at the upper an
lower bands respectively.

The residual contribution of the ionosphere is caused by a) higher order terms
in the expansion of the dispersiveness on frequency (Hawarey et al., 2005); b)
the contribution of frequency-dependent source structure, and c) the
dispersiveness in the signal chain. These contributions affect group delay at
a level of several picoseconds and they are considered insignificant with
respect to other systematic errors.

###  3.2 The use of TEC maps for computation of ionospheric path delay

TEC maps (Hernández-Pajares et al., 2009) also known as Global Ionospheric
Model (GIM) derived from analysis of Global Navigation Satellite System (GNSS)
data are used for reduction of single band observations. In particular, I used
CODE TEC time series (Schaer, 1999) available at ftp://ftp.aiub.unibe.ch/CODE
since January 01 1995 with a spatial resolution of $5^{\circ}\times
2.5^{\circ}$. Time resolution varied. It was $24^{h}$ since 01 January 1995
through March 27, 1998; $2^{h}$ since March 28, 1998 through October 18, 2014;
and $1^{h}$ after that date. The ionosphere is considered as a thin shell at
the constant height $H_{i}$ of 450 km above the mean Earth’s radius. The
ionospheric contribution is expressed via TEC as

$\displaystyle\uptau_{i}=\frac{\displaystyle\alpha}{\displaystyle
f^{2}_{e}}\>M(e)\>\mbox{TEC},$ (6)

where $M(e)$ is the so-called thin shell ionospheric mapping function

$\displaystyle M(e)=\frac{\displaystyle
1}{\displaystyle\sqrt{1-\left(\frac{\displaystyle\bar{R}_{\oplus}}{\displaystyle
R_{\oplus}+H_{i}}\right)^{2}\,\cos^{2}{e_{\rm gc}}}},$ (7)

and $R_{\oplus}$ is the mean Earth’s radius and $e_{\rm gc}$ is the geocentric
elevation angle with respect to the radius vector between the geocenter and
the station. Computation of TEC value at a given moment of time is reduced to
computation of the position of the ionosphere piercing point at a given
azimuth and elevation and interpolation of TEC at the piercing point at that
latitude, longitude, and time.

TEC maps from GNSS is a coarse model of the ionosphere. Errors of $\uptau_{i}$
computed according to expression 5 are greater than the residual ionosphere
contribution of ionosphere-free linear combinations of group delays.
Therefore, dual-band delay observables are preferred when they are available.
However, there are two cases when they are not available: a) dual-band
observing sessions with some source detected only at one band; b) single-band
observing sessions. In these two cases we resort to computation of the
ionospheric contribution to path delays using GNSS TEC maps and evaluation of
uncertainties of these contributions.

###  3.3 Ionospheric contribution in dual-band observing sessions when a
source is detected at one band only

The simplest way to deal with a mixture of dual- and single- band data is to
process experiments three times: 1) using dual-band data of those observations
that detected a source in both bands, 2) using a low band data, and 3) using
an upper band data with applying ionospheric path delay computed from GNSS TEC
maps. However, typically only a fraction, 2 to 20% of observations is detected
at only one band; the rest of observations are detected at both bands.
Therefore, we can use available dual-band observations at a given observing
session to improve the TEC model.

I represent ionospheric path delay at stations $j,k$ as

$\displaystyle\begin{array}[]{lcl}\uptau_{i}(t)&=&b_{j}(t)-b_{k}(t)\>+\\\
&&\frac{\displaystyle\alpha}{\displaystyle f^{2}_{e}}\biggl{(}\Bigl{(}\rm
TEC_{j}(\phi_{j},\lambda_{j},t)+a_{j}(t)\Bigr{)}M(e_{j})\>-\\\
&&\phantom{aaaa}\Bigl{(}\rm
TEC_{k}(\phi_{k},\lambda_{k},t)+a_{k}(t)\Bigr{)}M(e_{k})\biggr{)},\end{array}$
(11)

where $b_{j}(t)=\sum_{i}^{n}B^{0}_{i}(t)b_{ij}$ is a delay bias expanded over
the B-spline basis of the 0th degree,
$a_{j}(t)=\sum_{i}^{n}B^{3}_{i}(t)a_{ij}$ is the TEC bias expanded over the
B-spline basis of the 3rd degree. $\phi,\lambda$ are coordinates of the
ionosphere piercing point that depend on positions of observing stations as
well as azimuths and elevations of observed sources.

The clock bias occurs due to path delay in VLBI hardware that is different at
different bands. This bias is constant for most of the experiments, however
occasionally breaks may happen at some stations. Epochs of these breaks
coincide with the epochs of breaks in the clock function. Expansion over the
B-spline basis of the 0th degree accounts for these breaks. (B-spline of the
0th degree is 1 within the range of knots [i, i+1] and 0 otherwise.)

I estimated parameters $a_{j}$ for all the stations and $b_{j}$ for all the
stations except the one taken as a reference using all available dual-band
observations of a given experiment using least squares with weights

$\displaystyle w_{i}=\frac{\displaystyle
1}{\displaystyle\sqrt{y^{2}+\frac{f_{u}^{4}\sigma^{2}(\uptau_{u})+f_{l}^{4}\sigma^{2}(\uptau_{l})}{(f_{u}^{2}-f_{l})^{2}}}},$
(12)

where $\sigma(\uptau_{u})$ and $\sigma(\uptau_{u})$ are group delay
uncertainties and $y$ is the error floor, 12 ps, introduced to avoid
observations with very high signal to noise ratios to dominate the solution.

The time span of the B-spline knot sequence for TEC bias in my solutions was
15 minutes. I applied constraints on the value of the B-spline coefficients
and on the first and the second derivative with reciprocal weights $5\cdot
10^{-10}$ $s$, $4\cdot 10^{-14}$, and $2\cdot 10^{-18}$ $s^{-1}$ respectively.
These constraints were applied to ensure the continuity of biases and to
prevent a singularity in rare cases when too few available observations at a
given station could be used for bias estimation at a given spline segment.
Figure 1 illustrates estimates of the ionospheric bias.

Figure 1: Adjustment to the ionosphere path delay bias at 8.4 GHz with respect
to the path delay derived from GNSS TEC maps at mk-vlba station from
processing dual-band observations on April 22, 2015. 1 TECU causes group delay
21 ps at 8.4 GHz.

The resulting total electron contents model $\rm TEC_{j}(t)+a_{j}(t)$ is more
precise than the a priori $\rm TEC_{j}(t)$ taken from GNSS maps because it
uses additional information. Using estimates of $a_{j}$ and $b_{j}$ spline
coefficients, I compute $\uptau_{i}(t)$ and its uncertainty according to the
law of error propagation using the full variance-covariance matrix of spline
estimate coefficients. In order to evaluate the realism of these errors, I
processed the trial dataset and computed $\uptau_{i}(t)$ using the estimates
of clock and TEC biases and compared them with the ionospheric contribution
derived from dual-band observations. I removed clock biases from VLBI dual-
band ionospheric contributions $\uptau_{\rm vi}$, formed the differences
$\uptau_{i}-\uptau_{\rm vi}$, and then divided them by $\sigma(\uptau_{i})$
derived from the variance-covariance matrix of $a_{j},b_{j}$. I generated the
normalized histogram from the dataset of 4,343,782 differences and computed
the first two moments of the empirical distribution shown in Figure 2. The
fitting parameters of the first and second moments of the distribution are
0.003 and 0.889 respectively. Two factors cause a deviation of the second
moment from 1.0 in the opposite way: a) TEC variations not accounted by the
parametric model; b) statistical dependence of the estimates of $a_{j},b_{j}$
and VLBI path delay. After scaling the variance-covariance matrix by square of
$0.889$, the distribution of the normalized residuals becomes close to
Gaussian. The closeness of the empirical distribution to the normal
distribution provides us a confidence that the extra noise introduced by the
mismodeled ionospheric path delay after applying clock and TEC biases is
properly accounted for.

Figure 2: Empirical distribution of the normalized differences of the
ionosphere path delay computed from the GNSS TEC maps adjusted for clock and
TEC biases (green dots). The normal distribution with $\sigma=1$ is shown as a
reference (solid blue line).

The closeness of the distribution of normalized differences to Gaussian with
$\sigma=1$ is encouraging, but it does not guarantee that the residual errors
of the sum of TEC from GNSS and TEC bias adjustments cause no systematic error
in estimates of source position. To characterize the impact of residual errors
of the ionospheric contribution on source position, I ran a solution XIA that
had the following differences with respect to the reference solution: 1) it
used X-band group delays; 2) data reduction for the ionosphere accounted for
both a priori TEC from GNSS maps and the ionosphere bias adjustment
(expression in the parentheses 11); 3) errors of the ionospheric biases
$\sigma=a_{\rm ij}\,\alpha/f^{2}_{e}$ were added in quadrature to reciprocal
weights of observables. I ran also solution SIA that differed from XIA by
using S-band group delays and the S-band effective ionospheric frequency
$f_{e}$. For control, I ran solutions XIN and SIN that used the ionosphere-
free combinations of group delays, the same weights as XIA and SIA solutions
respectively, and the zero mean Gaussian noise with $\sigma=a_{\rm
ij}\,\alpha/f^{2}_{e}$ was added to each observable. The differences in XIA
and SIA solutions with respect to the reference solution provide us a measure
of the impact of residual ionospheric errors on source position. The
difference in XIN and SIN solutions characterize the impact of the residual
ionospheric errors if they were Gaussian and totally uncorrelated.

Figure 3 shows the distribution of differences in declinations of source
position estimates from XIA and XIN solutions. There is no noticeable
deviation from the Gaussian shape. Table 1 lists first two moments of the
distribution of position differences in both right ascension and declination.
The second moment from SIN solution is close to 1.0, while the second moment
from XIN solution is 0.56. The errors of the ionospheric contribution at
S-band dominate the error budget. These errors are 14 times less at X-band and
are only a fraction of overall group delay errors. The mean biases in right
ascension and declination are negligible. The second moment of position
estimates from XIA and SIA solutions is 15% higher than moments from positions
from XIN and SIN solutions respectively. This increase occurs due to non-
randomness of residual ionospheric errors and can be viewed as a measure of
unaccounted systematic errors in source positions due to the ionosphere.
Analysis of these trial solutions demonstrates that we are able to predict the
impact of ionospheric errors on source position with an accuracy of 15% when
TEC bias is estimated using dual-band observations.

Parameter | TEC biases adjusted | Gaussian noise added
---|---|---
| mean | $\sigma$ | mean | $\sigma$
$\Delta\alpha$ X-band | -0.03 | 0.63 | -0.01 | 0.56
$\Delta\delta$ X-band | 0.02 | 0.64 | 0.01 | 0.56
$\Delta\alpha$ S-band | -0.05 | 1.14 | -0.01 | 1.03
$\Delta\delta$ S-band | 0.06 | 1.13 | 0.02 | 1.02
Table 1: The first and second moments of position differences of trial
solutions XIA, SIA, XIN, and SIN with respect to the position estimates
derived from the reference solution. The a priori TEC from GNSS and bias
adjustment from dual-band solutions were applied in XIA and SIA solutions (the
2nd column). The zero-mean Gaussian noise with $\sigma$ equal to the
uncertainty in path delay from the TEC bias adjustment was added in XIN and
SIN solutions (the 3rd column).

Figure 3: The distribution of normalized differences in declination from the
trial VLBI solutions with respect to the reference dual-band solution. Left:
solution with GNSS TEC maps + adjustments of TEC biases using dual-band VLBI
group delays (XIA). Right: solution with added Gaussian noise with $\sigma$
equal to errors in path delay that correspond to TEC bias adjustments
uncertainties (XIN). Blue thick lines show fitted Gaussian distributions with
the first and second moments listed in Table 1.

###  3.4 Ionospheric contribution in single-band observing sessions

When an entire session is observed only at one band, TEC biases cannot be
computed. Therefore, we have to resort to deriving a regression model to
provide estimates of these errors. In the past, Petrov et al. (2019) derived a
regression against the so-called global TEC: the integral of TEC over the
entire Earth following ideas of Afraimovich et al. (2008) and Krásná & Petrov
(2021) derived a regression against the root mean square (rms) of the total
ionospheric path delay from GNSS TEC. In this study I use the second approach
with slight modifications. Following general results of the turbulence theory
(see Tatarskii, 1971), we can expect that fluctuations at scales $x$ are
related to fluctuations at scales $y$ via a power law.

I processed the same dataset of 263 twenty-four hour VLBI experiments that was
used in the previous section and computed residual ionospheric path delay for
each observation as

$\displaystyle\uptau_{r}=(\uptau_{\rm gi}-\uptau_{\rm
vi}-(b_{j}-b_{k}))\;\tilde{M},$ (13)

where $\uptau_{\rm gi}$ is the vertical ionospheric path delay from GNSS TEC
maps, $\uptau_{\rm vi}$ is the vertical ionospheric path delay from VLBI,
$b_{j}-b_{k}$ are contributions of the clock bias, and $\tilde{M}$ is the
averaged ionosphere mapping function between stations 1 and 2 of a given
baseline: $\tilde{M}=(M(e_{1})+M(e_{2}))/2.0$. The clock biases are routinely
adjusted during analysis of VLBI observations and therefore, their
contribution on VLBI results, such as source positions, is entirely
eliminated. Subtracting them in expression 13, I eliminate their impact on
statistics as well. I used only twenty-four hour VLBI experiments for deriving
statistics because the ionospheric path delay strongly depends on Solar time,
especially at low latitudes, and statistics derived at shorter time intervals
are not representative.

Table 2: Coefficients of the B-spline expansion of the dependence of the rms of residual ionospheric path delay derived from GNSS TEC maps on the rms of the total ionospheric path delay from GNSS TEC maps at 8 GHz.  knot index | knot argument | B-spline value
---|---|---
| (ps) | (ps)
-2 |  | 6.3
-1 |  | 14.8
-0 |  | 23.5
-1 | 0.0 | 114.0
-2 | 35.0 | 114.0
-3 | 120.0 | 114.0
-4 | 1300.0 |

Figure 4 shows the dependence of rms of residuals $\uptau_{r}$ on the rms of
total ionospheric path delay from GNSS TEC maps $\uptau_{\rm gi}$. Each point
on the plot corresponds to the rms for a given baseline and a given observing
session. I confirm the early result of Krásná & Petrov (2021) but here I used
a much larger dataset. The result reported in Krásná & Petrov (2021) was
slightly affected by an error in computation of ionospheric group delays for a
case when some data are flagged because of radio interference. This error has
been fixed and the affected experiments have been reprocessed from scratches.
This dependence can be coarsely described as a square root of the rms of the
total ionospheric path delay. For a better approximation I sought a regression
in the form of an expansion over B-splines of the 3rd degree. The spline
coefficients computed using least squares are listed in Table 2.

Figure 4: Dependence of the rms of residual ionospheric path delay derived
from GNSS TEC maps $\sigma_{\rm gr}$ on the rms of the total ionospheric path
delay from these maps $\sigma_{\rm gt}$. No adjustment to TEC has been
applied. Path delay is computed for the reference frequency 8 GHz. The blue
smooth line shows the regression model in a form of a B-spline that fits the
data.

Using that regression, I developed the following algorithm for computation of
errors of the ionospheric path delay from GNSS TEC. First, coordinates of $K$
points uniformly distributed over the sphere are computed using a random
number generator. Then for each baseline and each time epoch azimuth and
elevation angles of that point are computed at both stations of the baseline,
and if the elevations above the horizon are greater than $5^{\circ}$ at both
stations, that point is selected for further computations. If not, the next
point is drawn. Then total ionospheric path delay
$\uptau_{i}(A_{1},e_{1},A_{2},e_{2})$ is computed using GNSS TEC maps. It is
worth mentioning here that unlike to troposphere path delay,
$\uptau_{i}(A,e)\neq\uptau_{i}(A,\pi/2)\cdot M(e)$, since path delay depends
on positions of the ionosphere piercing points. It is not sufficient to
compute the ionospheric path delay in zenith direction and then map it via
$M(e)$: latitude and longitude of the piercing point can be as far as 1000 km
from the station. Following this approach, we sample piercing points uniformly
distributed within the mutual visibility zone. The process is repeated for
1440 time epochs that cover the time interval of VLBI experiment under
consideration with a step of 1 minute. Then for each baseline
$\sigma(\uptau_{\rm gt})$ is computed over a time series of 1440 $\uptau_{i}$
values. Finally, the estimate of the rms of residual ionospheric path delay
derived from GNSS TEC maps is computed from this regression via rms of the
total ionospheric path delay as

$\displaystyle\sigma_{\rm rr}=\sum_{k}^{n}B^{3}_{k}(\sigma(\uptau_{\rm
gt}))\,\sqrt{M^{2}(e_{1})+M^{2}(e_{2})}.$ (14)

Baseline-dependent datasets are considered independent for this computation:
the mutual visibility at all the stations of the network at a given moment of
time is not enforced. For several baselines longer than 96% Earth diameter,
this algorithm has a poor performance for selecting points above $5^{\circ}$.
Therefore a minor modification is made for such an extreme case: the elevation
angle is fixed to $5^{\circ}$, mutual visibility is not enforced, and azimuths
are selected randomly within a range of $[0,2\pi]$ independently for both
stations.

In order to evaluate the validity of this regression model of residual
ionospheric path delay errors, I computed $\uptau_{r}$ from dual-band
observations, $\sigma_{\rm gt}$ following the algorithm described above for
263 twenty-four experiments, and then computed the histograms of normalized
residuals $\uptau_{r}/\sigma_{\rm rr}$. The histogram is presented in Figure
5. The first two moments of the distribution are -0.083 and 1.214
respectively. Since regression $\sigma_{\rm rr}$ was found using least
squares, the number of observations with $\sigma(\uptau_{\rm gr})$ less and
greater than $\sigma_{\rm rr}$ for given $\sigma(\uptau_{\rm gt})$ is
approximately equal: the thick blue line cuts the cloud of green points in
Figure 5 almost by half. However, the variance of the contribution of those
points with $\sigma(\uptau_{r})>\sigma_{\rm rr}$ overweights the contribution
of those points with $\sigma(\uptau_{r})<\sigma_{\rm rr}$ because variance
quadratically depends on $\uptau_{r}$. This causes a positive bias. After
multiplying $\sigma_{\rm rr}$ by 1.214, the distribution of normalized
residuals becomes almost Gaussian.

One can notice that $\sqrt{M^{2}(e_{1})+M^{2}(e_{2})}$ in expression 14 is not
the same as $\tilde{M}=(M(e_{1})+M(e_{2}))/2$ used for computation of the
regression. I found that using $\tilde{M}$ instead of
$\sqrt{M^{2}(e_{1})+M^{2}(e_{2})}$ decreases the second moment from to 1.214
to 1.196, which is negligible,

Figure 5: The distribution of the normalized differences of the ionospheric
path delay computed from the GNSS TEC maps against VLBI ionospheric path delay
with clock biases subtracted (green dots). The normal distribution with
$\sigma=1$ (solid blue line) is shown as a reference.

The distributions shown in Figures 2 and 5 are computed for the entire dataset
of 4.3 million path delays and they represent the general population over the
interval of 23 years. Statistics for an individual observing session may
differ. In order to evaluate the scatter of the statistics, I computed the
time series of second moments of the distribution of normalized residuals of
ionospheric path delays and their uncertainties with and without TEC biases
adjusted for each observing session separately. I divided the normalized
residuals by scaling factors of 0.889 and 1.196 respectively. I computed the
distribution of second moment estimates and showed it in Figure 6. The scatter
of the second moments is small when TEC biases are adjusted. That means this
statistics is robust. When TEC is not adjusted, the scatter is significantly
larger, but even in that case 90% of the second moment estimates deviate from
1.0 by no more than 30%. This provides us a measure of uncertainties in
computation of ionospheric path delay errors in individual single-band
observing sessions.

Figure 6: The distribution of second moment estimates of the normalized
differences of ionospheric path delays derived from VLBI dual-band
observations and GNSS TEC maps among individual observing sessions. The narrow
green curve shows the statistics of the normalized residuals with TEC biases
adjusted and the wide blue curve show the statistics of normalized residuals
without TEC adjustment.

###  3.5 The impact of the residual ionospheric errors in source position in a
case of single-band observing sessions

Error analysis presented in the previous section characterizes our ability to
predict the first and second moments of the distribution, but it does not
guarantee that residual errors due to the ionosphere do not cause systematic
errors in source positions. I ran trial solutions XIT and SIT that used X- and
S-band group delays respectively, applied ionospheric path delays derived from
GNSS TEC maps, and inflated reciprocal weights of observables by adding in
quadrature errors in ionospheric path delays from TEC maps.

Analysis of differences in source positions from trial XIT and SIT solutions
with respect to the reference dual-band solution revealed no peculiarities in
right ascensions, but revealed significant systematic errors in declinations
(see Figure 7). The pattern of systematic errors in declination from the
S-band solution is similar but greater by a factor of
$f^{2}_{x}/f^{2}_{s}\approx 14$. We cannot consider a solution with such
errors as satisfactory. This was unexpected because prior work of Sekido et
al. (2003); Hobiger et al. (2006); Dettmering et al. (2011); Motlaghzadeh et
al. (2022) claimed a good agreement between TEC derived from VLBI and GNSS.
And indeed, the plot of ionospheric contributions in zenith direction from
VLBI after removal clock biases against the ionospheric contributions from
GNSS (Figure 8) shows no peculiarities and fits the straight line $\uptau_{\rm
vi}=-3.9\>{\rm ps}+1.06\cdot\uptau_{\rm gi}$. Although the residuals of this
dependence look random, they still cause systematic errors in source position.

Figure 7: Differences in declinations from the X-band single band solution XIT
with the data reduction for the ionosphere using ionospheric path delay from
GNSS TEC maps applied with respect to the dual-band reference solution. The
thick blue line shows the differences smoothed with the Gaussian filter.
Figure 8: Dependence of VLBI ionospheric group delay at 8 GHz against the
ionospheric group delay from GNSS. The blue straight line is the least square
fit of this dependence.

I made a number of trial solutions. The following leads turned out productive:
1) to modify mapping function and 2) to scale ionospheric path delays from
TEC. Schaer (1999) suggested the following modification of the ionospheric
mapping function:

$\displaystyle M(e)=k\,\frac{\displaystyle
1}{\displaystyle\sqrt{1-\left(\frac{\displaystyle\bar{R}_{\oplus}}{\displaystyle
R_{\oplus}+H_{i}+\Delta H}\right)^{2}\,\cos^{2}{\alpha e_{\rm gc}}}},$ (15)

arguing that varying parameters $\alpha$ and $H_{i}$ one can account for a
more realistic electron density distribution with height than a thin shell
model. Here $\Delta H$ is an increment in the ionosphere height and $\alpha$
is a fudge factor.

I ran trial solutions: XIE0 with $k=1.0,\,\Delta H=0,\,\alpha=1.00$, XIE1 with
$k=1.0,\,\Delta H=56.7$ km, $\alpha$ = 0.9782, and XIE2 with $k=1.0,\,\Delta
H=150.0$ km, $\alpha$ = 0.9782. The choice of $k,\Delta H,\alpha$ in XIE0
solution corresponds to the thin shell model and the choice in XIE1 solution
corresponds to the so-called JPL modified single layer ionospheric mapping
function. That flavor of the mapping function was used in a number of papers
(Li et al., 2019; Xiang & Gao, 2019; Zhao et al., 2021; Wielgosz et al.,
2021), however no details how these values of $k,\,\Delta H,$ and $\alpha$
were derived have been provided.

I estimated admittance factors of the ionospheric path delay from GNSS maps
using all X-band group delays in a number of least square solutions. No
ionospheric contribution was applied in these trials. The partial derivative
of the admittance factor over observed group delay is the GNSS ionospheric
path delay itself. Table 3 shows the estimates of the admittance factors.

Table 3: Estimates of the admittance factors of the ionospheric path delay from GNSS TEC maps from X-band group delays.  Solution | Admittance
---|---
XIE0 | $0.750\pm 0.002$
XIE1 | $0.802\pm 0.002$
XIE2 | $0.835\pm 0.002$

The admittance factor estimates deviate from 1.0 at a level of
$\sim\\!\\!100\sigma$, i.e. on an extremely high significance level. The
admittance factor becomes closer to 1.0 when parameter height of the
ionosphere in the mapping function is increased. I interpret this result as a
deficiency of GNSS TEC maps and I associate the origin of this deficiency with
an over-simplification of the mapping function that was used for derivation of
the TEC maps from processing of GNSS observations.

Unfortunately, it does not look possible to fix the deficiency of the GNSS TEC
maps without re-processing GNSS observations, which is well beyond the scope
of this article, but nevertheless, it is still feasible to mitigate the impact
of the imperfection of GNSS TEC maps on source position estimates.

I ran 12 trial solutions with mapping functions with parameters 1) $\Delta
H=0.0,\,\alpha=1.0$; 2) $\Delta H=56.7~{}{\rm km},\,\alpha=0.9782$; and 3)
$\Delta H=150.0~{}{\rm km},\,\alpha=0.9782$. I varied the mapping function
scaling factor $k$ setting it to 0.7, 0.8, 0.9, and 1.0 and computed
declination biases with respect to the reference solution. The results are
presented in Figures 9–11.

Figure 9: Smoothed declination biases from four solutions with respect to the
dual-band reference solution using X-band observables and the ionospheric
mapping function with parameters $\Delta H=0$, $\alpha=1.0$ with different
scaling factors $k=0.7,\,0.8,\,0.9,\,{\rm and}1.0$.  Figure 10: Smoothed
declination biases from four solution with respect to the dual-band reference
solution using X-band observables and the ionospheric mapping function with
parameters $\Delta H=56.7$ km and $\alpha=0.9782$ with different scaling
factors $k=0.7,\,0.8,\,0.9,\,{\rm and}1.0$.  Figure 11: Smoothed declination
biases from four solutions with respect to the dual-band reference solution
using X-band observables and the ionospheric mapping function with parameters
$\Delta H=150$ km and $\alpha=0.9782$ with different scaling factors
$k=0.7,\,0.8,\,0.9,\,{\rm and}1.0$.

Close analysis of these figures reveals that the declination bias has three
components: 1) a constant bias; 2) a linear increase in the declination bias
with a decrease of declination; 3) a feature
$\delta(x\,\sin\delta+y\,\cos\delta)$ with the minimum at approximately
$-12^{\circ}$, where $\delta$ is declination. All three components depend on
parameters of the mapping function $\Delta H,\,\alpha$ and on scaling factor
$k$. Selecting appropriate $\Delta H,\,\alpha,\,k$, we can reduce declination
biases. I selected $\Delta H=56.7~{}{\rm km},\,\alpha=0.9782,$ and $k=0.85$.
This choice makes the weighted mean declination bias over all sources 0.013
mas. There is an element of subjectivity in this specific choice, since there
exist another combinations of $\Delta H,\alpha,k$ that provide the zero mean
bias. Wielgosz et al. (2021) showed that the mean rms of the difference
between GNSS TEC maps and GNSS slant TEC was reduced from 2 to 6% when this
modified mapping function was used. That choice of $\Delta H$ and $\alpha$ led
me to a selection of the specific scaling factor $k$.

Although scaling and modification of the ionospheric mapping function results
in a substantial reduction of the declination bias, the remaining bias is
still worrisome. To mitigate the bias even further, I introduce the ad hoc
correction for the ionospheric bias in the data reduction model. I smoothed
the biases $D(\delta)$ with the Gaussian filter with $\sigma=8^{\circ}$ —
these smoothed biases are shown in Figures 9–11 — and expanded them over the
basis of B-splines of the 3rd degree with 13 equi-distant knots in the range
of $[-50^{\circ},90^{\circ}]$. I added correction

$\displaystyle\frac{\partial\uptau}{\partial\delta}\,D(\delta)\,f^{2}_{d}/f^{2}_{e}$
(16)

to observables. Here $f_{d}$ is the reference frequency of the bias model and
$f_{e}$ is the effective ionospheric path delay of a given observation. Figure
12 shows the effect of applying de-bias correction. The bias has gone.

I ran two solutions using X-band only data (XIS) and S-band only data (SIS),
applied ionospheric path delay derived from the a priori TEC from GNSS maps
using the modified mapping function with parameters $\Delta H=56.7{\rm
km},\,\alpha=0.9782,\,k=0.85$, and applied the declination bias correction.
The reciprocal weights of observables were adjusted by adding in quadrature
the errors of residual ionospheric errors $\sigma_{\rm rr}$ modeled according
to regression expression 14. The first and the second moments of the
normalized differences in sources positions are presented in Table 4. The
second moments are less than 1.0 in the X-band solution. This indicates that
the ionospheric contribution is not the dominating error source.

Parameter | X-band solution | S-band solution
---|---|---
| mean | $\sigma$ | mean | $\sigma$
$\Delta\alpha$ | -0.07 | 0.45 | -0.11 | 0.82
$\Delta\delta$ | -0.02 | 0.50 | -0.06 | 1.00
Table 4: The first and second moments of position differences of trial
solutions XIS and SIS with respect to the position estimates derived from the
reference solution. The a priori ionospheric contribution was computed from
GNSS TEC maps using the modified mapping function with parameters $\Delta
H=56.7{\rm km},\,\alpha=0.9782,\,k=0.85$, the declination bias correction was
applied, but no reduction for TEC bias adjustment has been applied.

Figure 12: The declination bias with the ionospheric path delay using mapping
function with parameters $\Delta H=56.7{\rm km},\alpha=0.9782,k=0.85$ with
(Right) and without (Left) applying the empirical de-bias correction.

##  4 Discussion

###  4.1 Inaccuracy of GNSS TEC maps

A number of authors compared the ionospheric contribution from GNSS, VLBI,
DORIS, radio occultation observation, and dual-band satellite altimeters
(Sekido et al., 2003; Hobiger et al., 2006; Dettmering et al., 2011;
Hernández-Pajares et al., 2017; Cokrlic et al., 2018; Li et al., 2018, 2019;
Xiang & Gao, 2019; Wielgosz et al., 2021; Zhao et al., 2021; Motlaghzadeh et
al., 2022). A message these publications convey is there is a reasonable
agreement between GNSS TEC maps and other techniques and there are no majors
problems. Therefore, large systematic errors driven by mismodeling of the
ionospheric contribution derived from GNSS maps came as a big surprise. Do
results presented in this study contradict to prior publications?

Comparison in the past was often performed using very short continuous dataset
of VLBI observations, from 5 to 15 days (Hobiger et al., 2005; Dettmering et
al., 2011; Etemadfard et al., 2021; Motlaghzadeh et al., 2022). The level of
the agreement was characterized in terms of additive errors of GNSS TEC model.
Unfortunately, when data from a short time interval are analyzed, a
distinction between additive and multiplicative errors becomes blurry.
Characterizing the differences in terms of biases and rms of their scatter did
not turn productive in revealing of systematic errors. Moreover, an
unconscious bias of focusing a study on assessment of an agreement rather than
investigation of disagreements, that were just noticed and reported, diverted
attention from investigation the differences in depth.

However, reading between lines of published papers, we can find pieces of
evidences supporting findings in this work. The distribution of differences in
VLBI vertical TEC with respect TEC derived from GNSS TEC maps presented in
Figure 6 in Hobiger et al. (2006) shows a very significant skew. The
distribution has a negative mean and a much greater left tail than the right
tail. VLBI vertical TEC from that study appeared less than the TEC from GNSS
TEC maps, in agreement with what I have found. The authors did not investigate
that stark deviation of the distribution of over one million differences from
the Gaussian shape, only noting that the bias is less than 3 TECU. Considering
the errors are normally distributed and additive and considering that the
derivation of TEC from VLBI does not introduce serious errors, one can expect
to arrive to a Gaussian distribution of residuals. For instance, the
distribution of differences in Figure 5 indeed does not show any measurable
deviation from the Gaussian shape. However, if errors are multiplicative, the
residuals will be non-Gaussian and their distribution will be skewed.

Satellite altimetry using Jason satellites (Vaze et al., 2010, and references
therein) provides an independent way for assessment of a level of the
disagreement between direct vertical TEC measurements and GNSS TEC maps. Li et
al. (2018) showed that comparisons of the differences revealed significant
systematic biases that depend on geomagnetic latitude. However, an attempt to
characterize these additive biases in terms of rms was not very productive.
Liu et al. (2018) presented spatial distribution of the differences (Figures
6–7). That distribution strikingly reminds the average distribution of TEC
itself, suggesting the differences are multiplicative. A more recent paper of
Dettmering & Schwatke (2022) revealed a highly significant scaling factor
between GNSS TEC maps and four altimeter missions. The scaling factors varied
from 0.809 to 0.919 which is very close to what has been found in my analysis
of astrometric observations of extragalactic radio sources.

Li et al. (2019) performed a comparison of TEC maps with Jason satellite
altimetry and with ionospheric radio occultations (see the overview of this
technique in Bonafoni et al., 2019). They characterized the differences as a
superposition of an additive bias and a multiplicative scale factor. The scale
factors defined as TECJason/GIM and TECCOSMIC/TECGIM vary with time and
latitude and stay in a range of 0.7–0.9 at night and 0.9–1.0 during daytime.
Since Jason orbits have altitude 1,350 km and GNSS orbits have altitude about
20,000 km, Jason altimetry misses the contribution from the upper layers of
the ionosphere and plasmoshere. However, the electron density at altitudes
above 1,350 km is too low to account for 15–20% of the TEC that is not present
in Jason or COSMIC data. If the electron contents above altitude 1,350 km were
really so high to contribute to TEC at that level, the thin shell model with
the altitude of 450 km would be grossly inadequate.

Figure 13: The rms of the differences in right ascensions (upper row) and
declinations (lower row) when single band group delays are used with respect
to the positions from the dual-band reference solution. A priori ionospheric
path delay using GNSS TEC with the modified mapping function ($\Delta
H=56.7~{}{\rm km},\alpha=0.9782,k=0.85$) was used adjusting TEC biases from
dual-band observations (right) and without adjusting biases but with applying
the de-bias correction from expression 16 (left). The upper blue band shows
the differences for S-band, next red line shows the differences for C-band,
next green line shows the differences for X-band, and the bottom purple line
shows the differences for K-band.

It is essential to note that analysis of Jason and COSMIC data provides
estimates of vertical TEC, while analysis if GNSS observations provides slant
TEC that is converted to vertical TEC using the mapping function. The
dependence of the admittance factor estimate on the mapping function found in
this study suggests that a simple thin model may not be adequate. Schaer
(1999) discussed the dependence of the effective ionosphere height on solar
zenith angle. Xiang & Gao (2019) studied it in more details. The height of the
peak electron density has annual and diurnal variations. The latter variations
have an amplitude of about 100 km, being lower at daytime. They showed that
the instantaneous mapping function that accounts for the height of the
electron content maximum achieves 8% reduction of mapping errors. It should be
mentioned that if the effective height of the ionosphere is changed, the
latitude and longitude of the ionosphere piercing point for an observation
with a given elevation and azimuth will change as well, which will provide an
additional change in path delay.

These works strengthen argumentation in favor of that a thin shell model used
for derivation of TEC maps since mid 1990s is oversimplified. The effective
height of the ionosphere changes with time and latitude and therefore, a
realistic mapping function should also vary with time and latitude. Omission
of this complexity results in a deficiency of GNSS TEC maps that manifests in
multiplicative (scale) and additive (bias) errors that varies with time and
latitude. The non-linear dependence of the declination bias with declination
in a form of $\delta(x\,\sin\delta+y\,\cos\delta)$ may be caused by the
dependence of GNSS TEC scaling factor with latitude reported by Li et al.
(2019). Radio waves from southern hemisphere sources observed at the array
located in the Northern Hemisphere propagate through regions in the ionosphere
with systematically lower latitudes than from northern hemisphere sources.
Since a fixed mapping function is used for both computation of TEC maps and
computation of slant ionospheric path delay from these maps, no modification
of that fixed mapping function for data reduction of astronomical data is able
to account for this kind of complexity, but as it was shown earlier, it is
still possible to mitigate it. The remaining bias can be eliminated by
applying an empirical de-bias correction.

###  4.2 Omitted refinements of the ionosphere contribution

Etemadfard et al. (2021) processed 60 days of VLBI data and estimated not only
TEC for each sitet using B-splines of the 1st degree, but also TEC derivatives
over longitude and latitude that were considered constant over 24 hour
periods. They claim that estimating TEC partial derivatives decreases
discrepancies of TEC from VLBI to GNSS TEC maps by 36%. Inspired by this
results, I introduced estimation of latitude and longitude TEC gradients in a
form of B-spline in addition to TEC estimation using dual-band data. I
considered the weighted rms (wrms) of the differences between the
parameterized model of TEC adjustment to vertical TEC from dual-band data as a
metric of an improvement. I did not find a reduction of rms greater then
several percents and abandoned this approach. I should note this result does
not disprove findings of Etemadfard et al. (2021) since they used a different
metric.

Expression for the ionospheric contribution 1 is an approximation. Hawarey et
al. (2005) considered the impact of the higher order of expansion on group
delay, namely proportional to $f^{-3}$. They found that the maximum
contribution of the 3rd degree term varied from 3 to 9 ps at 8 GHz depending
on the baseline. This contribution is approximately one order of magnitude
less than the contribution of residual ionospheric errors after taking into
account GNSS TEC maps. It should be noted that the 3rd degree term affects
both single-band and ionosphere-free linear combinations of group delays at
two bands and therefore, cannot be retrieved in analysis of the differences
with respect to the dual-band reference solution.

###  4.3 The impact of remaining errors in modeling of the ionosphere on
source positions

Figure 14: An increase in source position errors in right ascension (Left) and
declination (Right) due to the residual ionospheric contribution at K-band (22
GHz) after applying data reduction based on GNSS TEC maps with the modified
mapping function ($\Delta H=56.7~{}{\rm km},\alpha=0.9782,k=0.85$)

As it was shown before, our ability to model the ionospheric contribution
using single-band observations is limited. After applying the de-bias data
reduction, the declination bias is virtually eliminated. The residual
ionospheric contribution causes additional random errors. In order to
investigate their impact, I ran four trial solutions. I used ionosphere-free
linear combinations of dual-band data and added to them the ionospheric
contribution derived from VLBI data scaled to the specific frequency of a
trial solution. Then I applied the data reduction for the ionosphere to this
dataset of modified observables using GNSS TEC maps as if I processed single-
band observations and estimated source positions. The differences in source
positions from these solutions with respect to the reference solutions are
interpreted as an impact of the residual ionospheric errors at different
frequencies. The declination dependence of the differences in right ascensions
and declination for four solutions for frequency 2.3 GHz (S-band), 4.3 GHz
(C-band), 8.4 GHz (X-band), and 23.7 GHz (K-band) is shown in Figure 13.

These plots help to quantify additional errors in source positions that would
arise if the dataset of 263 twenty-four hour experiments used for this study
would have been observed at one band only. We see that estimation of TEC
biases reduces the ionosphere-driven errors by a factor of 2–4. Still, even
after estimation of TEC biases, the declination dependence of the ionosphere-
driven errors on declination remains.

It is known that unaccounted source structure affects source position
estimates. In general, the contribution of source structure at higher
frequencies is less because jets become optically thin. Therefore, there was
an expectation that observing at high frequencies, such as 22–24 GHz (K-band)
one may obtain more precise source positions. So far, observational evidence
do not support that prediction (Lanyi et al., 2010; Charlot et al., 2020). A
detailed comparison of K-band absolute astrometry versus dual-band astrometry
of Karbon & Nothnagel (2019) did not reveal an improvement. Therefore, it is
instructive to see what is the contribution of the ionospheric errors at
K-band after applying data reduction from GNSS TEC maps. Figure 14 shows
additional errors in right ascensions and declinations due to residual
contribution of the ionosphere. We see that additional errors in right
ascension are about 0.1 mas. Declination errors of northern hemisphere sources
are at the same level. However, these errors grow with a decrease of
declination for sources in the southern hemisphere approximately linearly and
reach 0.3 mas at declinations $-40^{\circ}$. Therefore, the unmodeled
ionospheric contribution sets the error floor in position accuracy. This error
floor is not accounted in the ICRF3 catalogue, and the ionospheric
contribution is 300% greater than the noise floor at K-band adopted by Charlot
et al. (2020) according to their Table 6. The potential of K-band astrometry
cannot be utilized unless the accuracy of ionosphere modeling will be
substantially improved.

###  4.4 The rms of residual atmospheric errors

I investigated how the rms of the residual ionospheric contributions for 45
VLBA baselines for the dateset of 263 twenty-four observing sessions varied
with time. I computed three statistics for each observing session: 1)
$\uptau_{v}-b_{i}$, 2) $\uptau_{v}-\uptau_{g}-b_{i}$; and 3)
$\uptau_{v}-\uptau_{g}-a_{i}-b_{i}$, and re-scaled them to 8 GHz. Here
$\uptau_{v}$ and $\uptau_{g}$ are ionospheric path delays from VLBI and TEC
maps respectively, $a_{i}$ is the adjusted TEC bias, and $b_{i}$ is clock
bias. The statistics are shown in Figure 15. In order to improve readability,
the time series were smoothed using Gaussian kernel with parameter $\sigma=1$
year. The first statistics characterizes the impact of the total ionosphere on
group delay. The second statistics characterizes the rms of the impact of the
residual ionospheric errors on group delay after applying the a priori GNSS
TEC model. The final statistics characterizes the impact of residual errors
after applying the a priori GNSS TEC and adjusting TEC biases using dual-band
VLBI data.

Figure 15: The rms of the mean residual ionospheric contribution at VLBA
baselines at 8 GHz for three cases: 1) no a priori ionospheric contribution is
applied (upper red curve); 2) the a priori ionospheric path delay computed
using GNSS TEC maps is applied (middle blue curve); and 3) the a priori
ionospheric path delay computed using GNSS TEC maps are applied and TEC biases
are adjusted (lower green curve). 1 TECU causes group delay 21 ps at 8 GHz.

###  4.5 A combined use of dual-band and single-band data

The lack of biases and the realistic assessment of errors due to residual
ionospheric contribution allows us to use a mixture of dual-band and single-
band group delays in a single least square solution. Source positions derived
from single-band data are in general less precise than those derived from
dual-band data because an additional factor affects the uncertainties of group
delays. These additional errors can be computed and accounted for in deriving
weights of observables. An increase in uncertainties of group delay
observables propagates to uncertainties of source positions.

Combined analysis of a heterogeneous datasets provides a significant advantage
because all the data are used. We can fuse dual-band and single-band data in
one dataset and use it for estimation of source positions. A fused dataset
consists of observables of three types:

  1. 1.

Dual-band ionosphere-free linear combinations of group delays.

  2. 2.

Single-band delays of dual-band experiments detected at one band only. Data
reduction accounts for the ionospheric path delay computed from a sum of GNSS
TEC and TEC bias adjustments. The uncertainty of the bias adjustment is added
in quadrature to the group delay uncertainty for reciprocal weights of such
observables.

  3. 3.

Group delays from single-band experiments. Data reduction accounts for the
ionospheric path delay computed from GNSS TEC and applies the de-bias
correction. The uncertainty of the ionospheric model computed from the
regression model is added in quadrature to the group delay uncertainty for
reciprocal weights of such observables.

Solving for source coordinates using a fused dataset provides estimates with
minimum uncertainties and minimum correlations since it uses all available
data under condition that the frequency-dependent position biases due to
source structure are negligible. This condition may be violated for some
strong sources with significant structure. Analysis of position offsets
between single band and dual-band solutions will help to identify such sources
(see, for example, Petrov et al., 2011), but these cases are encountered
infrequently.

###  4.6 Future work

An improvement in modeling of ionospheric path delay when TEC biases are
adjusted is quite impressive. However, this requires using additional
information about the state of the ionosphere. Specifically, VLBI ionospheric
path delay from dual-band data was used for computation of these biases. This
information is not available in processing of single-band observations. Ros et
al. (2000) considered the use of dual-band GPS observations from collocated
receivers for analysis of radio astronomy observations. They have shown that
the GPS determination of TEC from ground receivers alone without the use of
GNSS TEC maps can be successfully applied to the astrometric analysis of VLBI
observations.

There are on-going efforts to install advanced GNSS receivers within a hundred
meters of each VLBA antenna. The use of geometry-free pseudo-ranges at 1.2 and
1.5 GHz in a similar way as I used VLBI ionosphere-free group delays for
adjustments of TEC biases promises a similar level of improvement.

##  5 Conclusions

Single-band group delays are affected by the contribution of the ionosphere.
This contribution is noticeable at frequencies below 30 GHz and becomes a
dominating source of errors at frequencies below 5–8 GHz. Compared with
astrometric solutions based on the use of the ionosphere-free linear
combination of dual-band observables, an absolute source position catalogue
based on a single-band solution is affected by additional random and
systematic errors caused by mismodeling the contribution of the ionosphere to
path delay. I explored two approaches to modeling the ionospheric path delay
using GNSS TEC and assessed residual ionospheric errors.

The findings can be summarized as follows:

  1. 1.

In a case when an experiment was recorded at two widely separated bands, but a
fraction of sources were detected at one band only, estimation of the TEC bias
in a form of an expansion over B-spline basis using dual-band data and then
applying that bias in addition to GNSS TEC provides an unbiased estimates of
sources positions. The stochastic model that describes residual errors of TEC
bias adjustment predicts an increase of positions errors with an accuracy of
15%. No remaining systematic errors were found. This approach provides source
positions with the lowest uncertainties with respect to other approaches.

  2. 2.

In a case of single-band observations, path delay is computed using GNSS TEC
maps. Using the thin spherical shell model at the ionosphere height 450 km
above the mean Earth’s radius causes a strong systematic bias in declination
that reaches 1 mas at 8 GHz an 12 mas at 2.3 GHz. This bias can be virtually
eliminated when a) the modified ionospheric mapping function with parameters
$\Delta H=56.7~{}{\rm km},\,\alpha=0.9782,\,k=0.85$ is used; and b) the
empirical de-bias correction is applied.

  3. 3.

Ionospheric errors are mainly multiplicative.

  4. 4.

The scaling factor of GNSS TEC that provides the zero mean average declination
bias, 0.85, is in a good agreement with a totally independent comparison of
vertical TEC determined from satellite altimetry with GNSS TEC maps. Since
VLBI is sensitive to a delay incurred in the total ionosphere, interpretation
of a scaling factor of TEC from altimetry and radio occultation as a
contribution of upper layers of the ionosphere at altitudes above Jason orbit,
i.e. 1,350 km, suggested by Liu et al. (2018) and Dettmering & Schwatke (2022)
is ruled out.

  5. 5.

I have found that the scaling factor of GNSS TEC maps that provides zero mean
declination bias depends on the used ionospheric mapping function. Therefore,
I surmise that the established deficiency of GNSS TEC maps is caused by over-
simplification of the ionospheric mapping function used for their derivation
that considers the ionosphere as a thin spherical shell. The electron contents
in the real ionosphere varies not only in latitude and longitude, but also in
height. Diurnal variations of the effective ionosphere height at a level of
100 km, i.e. over 20%, are large enough to cause errors of the magnitude that
was found in analysis of VLBI data.

  6. 6.

The impact of the ionosphere on path delay depends on the Solar cycle.
Modeling ionospheric path with GNSS TEC maps reduces the residuals at 8 GHz by
a factor of 2 during the Solar maximum and only by 10% during the solar
minimum. Estimation of the TEC bias reduces ionospheric errors further by a
factor of 2 regardless of the phase of the Solar maximum.

  7. 7.

The impact of the ionosphere on source position errors can be modeled with an
accuracy of 15%. It remains noticeable even at frequencies as high as 22–24
GHz (K-band). In particular, the ionospheric errors even after applying data
reduction based on GNSS TEC with the modified mapping function and the de-bias
correction exceed 0.1 mas. Declination errors of sources with declinations in
a range of [$-40^{\circ},0^{\circ}$] are in a range of 0.1 to 0.3 mas. An
assertion that K-band astrometry is able to provide results more precise than
0.1 mas is not true at the current state of our ability to model ionospheric
path delay. Considering on-going efforts to install advanced GNSS receivers in
the close vicinity of VLBA stations and other radio telescopes, the situation
may change in the future.

This study lays the foundation of the single-band absolute astrometry. Dual-
band astrometric observations still provide the best accuracy. The use of
singe-band data with the procedure of data reduction and weighting described
above allows us to get unbiased positions with known added errors.

This work was done with datasets RDV, RV, CN, UF001, UG003, BL122, BL166,
BP133, BP138, GC073, BC204, BG219, and V17 collected with VLBA instrument of
the NRAO and available at https://archive.nrao.edu/archive. The NRAO is a
facility of the National Science Foundation operated under cooperative
agreement by Associated Universities, Inc. This work made use of the Swinburne
University of Technology software correlator, developed as part of the
Australian Major National Research Facilities Programme and operated under
license. It is my pleasure to thank Yuri Y. Kovalev, Urs Hugentobler, and
Robert Heinkelmann for discussions of presented results. VLBA

## References

  * Afraimovich et al. (2008) Afraimovich, E. L., Astafyeva, E. I., Oinats, A. V., Yasukevich, Y. V., & Zhivetiev, I. V. 2008, Annales Geophysicae, 26, 335, doi: 10.5194/angeo-26-335-2008
  * Behrend et al. (2020) Behrend, D., Thomas, C., Gipson, J., Himwich, E., & Bail, K. L. 2020, Journal of Geodesy, 94, doi: 10.1007/s00190-020-01436-x
  * Bonafoni et al. (2019) Bonafoni, S., Biondi, R., Brenot, H., & Anthes, R. 2019, Atmospheric Research, 230, 104624, doi: 10.1016/j.atmosres.2019.104624
  * Charlot et al. (2020) Charlot, P., Jacobs, C. S., Gordon, D., et al. 2020, Astronomy & Astrophysics, 644, A159, doi: 10.1051/0004-6361/202038368
  * Cokrlic et al. (2018) Cokrlic, M., Wezka, K., Galas, R., et al. 2018, in 2018 International Symposium ELMAR (IEEE), doi: 10.23919/elmar.2018.8534598
  * Dettmering et al. (2011) Dettmering, D., Heinkelmann, R., & Schmidt, M. 2011, Journal of Geodesy, 85, 443, doi: 10.1007/s00190-011-0473-z
  * Dettmering & Schwatke (2022) Dettmering, D., & Schwatke, C. 2022, Earth and Space Science, 9, doi: 10.1029/2021ea002098
  * Etemadfard et al. (2021) Etemadfard, H., Heinkelmann, R., Hossainali, M. M., & Schuh, H. 2021, Space Weather, 19, doi: 10.1029/2020sw002633
  * Gordon et al. (2016) Gordon, D., Jacobs, C., Beasley, A., et al. 2016, AJ, 151, 154, doi: 10.3847/0004-6256/151/6/154
  * Hawarey et al. (2005) Hawarey, M., Hobiger, T., & Schuh, H. 2005, Geophys. Res. Lett., 32, L11304, doi: 10.1029/2005GL022729
  * Hernández-Pajares et al. (2017) Hernández-Pajares, M., Roma-Dollase, D., Krankowski, A., García-Rigo, A., & Orús-Pérez, R. 2017, Journal of Geodesy, 91, 1405, doi: 10.1007/s00190-017-1032-z
  * Hernández-Pajares et al. (2009) Hernández-Pajares, M., Juan, J. M., Sanz, J., et al. 2009, Journal of Geodesy, 83, 263, doi: 10.1007/s00190-008-0266-1
  * Hobiger et al. (2005) Hobiger, T., Boehm, J., & Schuh, H. 2005, in International Association of Geodesy Symposia, ed. F. Sanso, Vol. 128 (Springer-Verlag), 131–136, doi: 10.1007/3-540-27432-4_23
  * Hobiger et al. (2006) Hobiger, T., Kondo, T., & Schuh, H. 2006, Radio Science, 41, RS1006, doi: 10.1029/2005rs003297
  * Karbon & Nothnagel (2019) Karbon, M., & Nothnagel, A. 2019, Astronomy & Astrophysics, 630, A101, doi: 10.1051/0004-6361/201936083
  * Klobuchar (1996) Klobuchar, J. A. 1996, in Global Positioning System: Theory and Applications, Volume I, ed. J. J. S. Jr., P. Axelrad, B. W. Parkinson, & P. Enge (American Institute of Aeronautics and Astronautics), 485–515, doi: 10.2514/4.866388
  * Krásná & Petrov (2021) Krásná, H., & Petrov, L. 2021, Journal of Geodesy, 95, 101, doi: 10.1007/s00190-021-01551-3
  * Lanyi et al. (2010) Lanyi, G. E., Boboltz, D. A., Charlot, P., et al. 2010, AJ, 139, 1695, doi: 10.1088/0004-6256/139/5/1695
  * Li et al. (2018) Li, M., Yuan, Y., Wang, N., Li, Z., & Huo, X. 2018, GPS Solutions, 22, doi: 10.1007/s10291-018-0721-2
  * Li et al. (2019) Li, W., Huang, L., Zhang, S., & Chai, Y. 2019, Sensors, 19, 5489, doi: 10.3390/s19245489
  * Liu et al. (2018) Liu, L., Yao, Y., Kong, J., & Shan, L. 2018, Remote Sensing, 10, 621, doi: 10.3390/rs10040621
  * Motlaghzadeh et al. (2022) Motlaghzadeh, S., Alizadeh, M. M., Cappallo, R., Heinkelmann, R., & Schuh, H. 2022, Radio Science, 57, doi: 10.1029/2021rs007336
  * Petrov (2021) Petrov, L. 2021, AJ, 161, 14, doi: 10.3847/1538-3881/abc4e1
  * Petrov et al. (2019) Petrov, L., de Witt, A., Sadler, E. M., Phillips, C., & Horiuchi, S. 2019, MNRAS, 485, 88, doi: 10.1093/mnras/stz242
  * Petrov et al. (2009) Petrov, L., Gordon, D., Gipson, J., et al. 2009, Journal of Geodesy, 83, 859, doi: 10.1007/s00190-009-0304-7
  * Petrov et al. (2011) Petrov, L., Kovalev, Y. Y., Fomalont, E. B., & Gordon, D. 2011, AJ, 142, 35, doi: 10.1088/0004-6256/142/2/35
  * Reid & Honma (2014) Reid, M. J., & Honma, M. 2014, ARA&A, 52, 339, doi: 10.1146/annurev-astro-081913-040006
  * Rioja & Dodson (2020) Rioja, M. J., & Dodson, R. 2020, The Astronomy and Astrophysics Review, 28, doi: 10.1007/s00159-020-00126-z
  * Ros et al. (2000) Ros, E., Marcaide, J. M., Guirado, J. C., Sardón, E., & Shapiro, I. I. 2000, A&A, 356, 357.  https://arxiv.org/abs/astro-ph/0001503
  * Schaer (1999) Schaer, S. 1999, Geod.-Geophys. Arb. Schweiz, 59
  * Sekido et al. (2003) Sekido, M., Kondo, T., Kawai, E., & Imae, M. 2003, Radio Science, 38, n/a, doi: 10.1029/2000rs002620
  * Tatarskii (1971) Tatarskii, V. I. 1971, The effects of the turbulent atmosphere on wave propagation (U. S. Department Of Commerce)
  * Vaze et al. (2010) Vaze, P., Neeck, S., Bannoura, W., et al. 2010, in SPIE Proceedings, ed. R. Meynart, S. P. Neeck, & H. Shimoda (SPIE), doi: 10.1117/12.868543
  * Wielgosz et al. (2021) Wielgosz, P., Milanowska, B., Krypiak-Gregorczyk, A., & Jarmołowski, W. 2021, GPS Solutions, 25, doi: 10.1007/s10291-021-01142-x
  * Xiang & Gao (2019) Xiang, Y., & Gao, Y. 2019, Remote Sensing, 11, 1497, doi: 10.3390/rs11121497
  * Zhao et al. (2021) Zhao, J., Hernández-Pajares, M., Li, Z., Wang, N., & Yuan, H. 2021, Journal of Geodesy, 95, doi: 10.1007/s00190-021-01487-8


# KD-Lib: A PyTorch library for Knowledge Distillation, Pruning and
Quantization

Het Shah,1 Avishree Khare,2 Neelay Shah,3∗ Khizir Siddiqui 4∗ Equal
contribution

###### Abstract

In recent years, the growing size of neural networks has led to a vast amount
of research concerning compression techniques to mitigate the drawbacks of
such large sizes. Most of these research works can be categorized into three
broad families : Knowledge Distillation, Pruning, and Quantization. While
there has been steady research in this domain, adoption and commercial usage
of the proposed techniques has not quite progressed at the rate. We present
KD-Lib, an open-source PyTorch based library, which contains state-of-the-art
modular implementations of algorithms from the three families on top of
multiple abstraction layers. KD-Lib is model and algorithm-agnostic, with
extended support for hyperparameter tuning using Optuna and Tensorboard for
logging and monitoring. The library can be found at -
https://github.com/SforAiDl/KD˙Lib

## Introduction

Deep neural networks (DNNs) have gained widespread popularity in recent years,
finding use in several domains including computer vision, natural language
processing, human computer interaction and more. These networks have achieved
remarkable results on several tasks, often even surpassing human-level
performance.

The number of parameters of such DNNs often increase multi-fold with an
increase in their representation capacity, limiting the deployment
capabilities and hence, the commercial feasibility of these networks. This
limitation warrants the need for efficient compression techniques that can
shrink the networks in size while ensuring that the drop in performance is
minimal. In this paper, we restrict our focus to three widely-used compression
techniques: Knowledge Distillation, Network Pruning and Quantization.

Knowledge Distillation (Hinton, Vinyals, and Dean 2015) is a compression
paradigm that leverages the capability of large neural networks (called
teacher networks) to transfer knowledge to smaller networks (called student
networks). While large models (such as very deep neural networks or ensembles
of many models) have higher knowledge capacity than small models, this
capacity might not be fully utilized. It can be computationally just as
expensive to evaluate a model even if it utilizes little of its knowledge
capacity. Knowledge distillation aims to transfers knowledge from a large
model to a smaller model without loss of validity. Several advancements have
been witnessed in the development of richer knowledge distillation algorithms,
attempting to reduce the difference in test accuracies of the teacher and the
student. These algorithms are model-agnostic and hence can be used for a wide
variety of network architectures.

While knowledge distillation attempts to train an equally-competent smaller
network, network pruning (LeCun et al. 1990) attempts to reduce the size of
the existing network by removing unimportant weights. Different pruning
techniques differ in the choice of weights to eliminate and the methods used
to do the same. Pruning can help in reducing the size of the network up to 90%
with minimal loss in performance. Some approaches have also been empirically
shown to result in faster training of the pruned network along with a higher
test accuracy (Frankle and Carbin 2018).

Quantization is another way to compress neural networks by reducing the number
of bits used to store the weights. As the weights of a network are usually
stored as 32-bit floating values (FP32), reducing the precision to 8-bit
integer values (INT8) will reduce the size of the network by 4 times. Several
approaches have been developed to quantize networks with minimal loss in
performance.

These compression techniques have become extremely popular in recent years and
are actively being researched. New algorithms proposed in research papers can
be difficult to understand and implement, especially for potential users in a
non-academic setting, thereby limiting their commercial usage. To the best of
our knowledge, there does not exist an umbrella framework containing
implementations of state-of-the-art algorithms in Knowledge Distillation,
Pruning and Quantization. In this paper, we present KD-Lib, a comprehensive
PyTorch based library for model compression. KD-Lib aims to bridge the gap
between research and widespread use of model compression techniques. We
envision that such a framework would be helpful to researchers as well,
providing them a tool to build upon existing algorithms and helping them in
going from idea to implementation faster.

Library | Knowledge Distillation | Pruning | Quantization
---|---|---|---
KD-Lib (Ours) | Present | Present | Present
Distiller(Zmora et al. 2019) | Present (only 1 algorithm) | Present | Present
AIMET3 | - | Present | Present
AquVitae1 | Present | - | -
Distiller1 | Present | - | -
Table 1: Comparision of various libraries with KD-Lib

## Related work

We compare KD-Lib with several openly available frameworks and libraries. In
our comparison, we do not include libraries that support less than two
algorithms.

Distiller (Zmora et al. 2019) is the most extensive framework we found, but it
primarily focuses on quantization and pruning with only one knowledge
distillation algorithm (Hinton, Vinyals, and Dean 2015).
AquVitae111https://github.com/aquvitae/aquvitae contains 4 distillation
methods but no quantization and pruning algorithms. Similarly
Distiller222https://github.com/karanchahal/distiller has 11 knowledge
distillation techniques but lacks pruning and quantization methods.
AIMET333https://github.com/quic/aimet focuses mainly on quantization and some
other relatively less popular model compression techniques such as tensor
decomposition. In our survey, we found no library containing algorithms
pertaining to all 3 of the popular compression paradigms - knowledge
distillation, pruning and quantization. Table 1 shows concise comparison with
different frameworks.

## Features and Algorithms

KD-Lib houses several algorithms proposed in recent years for model
compression. The following features have driven the design choices for the
library:

  * •

The main aim of KD-Lib is to make model compression algorithms accessible to a
wide range of users, and hence the work is fully open-source.

  * •

The library should act as a catalyst for further research in these fields. It
should also be extendable to newer algorithms and other model compression
fields. Hence, it is designed to be modular, allowing flexible modifications
to essential components that can lead to novel algorithms or better extensions
to existing algorithms.

  * •

The interface should be easy to use. Hence, the core functionalities
(distillation/pruning/quantization) are accessible in a few lines of code.

  * •

As tuning the hyperparameters is essential for optimum performance, KD-Lib
provides support for hyperparameter tuning via Optuna. Monitoring and logging
support is also provided through Tensorboard.

A brief description of the implemented algorithms is as follows:

  * •

Knowledge Distillation : The algorithms have been divided into two major task-
types: Vision and Text. The Vision module currently supports 13 algorithms
while the Text module supports distillation from BERT to LSTM-based
networks(Tang et al. 2019).

  * •

Pruning : The library currently supports pruning based on the Lottery ticket
Hypothesis (Frankle and Carbin 2018).

  * •

Quantization : Static Quantization, Dynamic Quantization and Quantization
Aware Training (QAT) (Jacob et al. 2018) are currently supported by KD-Lib.

## Code Structure

The structure of the library has been designed for efficient use with the
following major principles kept in mind:

  * •

The core function of an algorithm can be executed in one line of code. Hence,
the classes contain a dedicated method for distillation/pruning/quantization.

  * •

Each module allows extension to newer features and easy modifications. Hence,
fluid components of algorithms (loss functions in distillation, for example)
can be easily customized.

  * •

Necessary statistics are available wherever needed. Hence, methods dedicated
to these are also present ($get\\_pruning\\_statics$, for example).

.1 Distiller. .2 train_student. .2 train_teacher. .2 evaluate. .2
calculate_kd_loss.

Figure 1: Structure of a Distiller.

Knowledge Distillation algorithms can be accessed as Distiller objects (Figure
1), with at least the mentioned methods. The train_student method distills
knowledge from a teacher network to a student network, where the teacher
network could optionally be trained using the train_teacher method. The
evaluate method can be invoked to test the performance of the student network.
The calculate_kd_loss method can overridden to provide a custom loss function
for distillation. This can also be leveraged by researchers to test novel
Knowledge Distillation loss functions.

.1 Pruner. .2 prune. .2 get_pruning_statistics.

Figure 2: Structure of a Pruner.

Pruning algorithms have been implemented as Pruner objects (Figure 2). Each
Pruner object can access the prune method for pruning the network.
Additionally, the get_pruning_statistics method can be used to obtain
information about the weights of the network after pruning (percentage of
network pruned, for example).

.1 Quantizer. .2 quantize. .2 get_performance_statistics. .2 get_model_sizes.

Figure 3: Structure of a Quantizer.

Quantization algorithms can be accessed via Quantizer objects (Figure 3). The
quantize method can be used for quantization (with differing implementations
for different algorithms). Additionally, the get_model_sizes method can be
used to compare sizes of the model before and after quantization and the
get_performance_statistics method can be used to compare test-times and error
metrics for the two networks.

The documentation for the library444https://kd-lib.readthedocs.io/ has the
description of all classes and selected tutorials with example code snippets.

## Benchmarks

We summarize benchmark results on some of the algorithms implemented in KD-Lib
in Tables 2, 3 and 4.

Algorithm | Accuracy
---|---
None | 0.57
DML (Zhang et al. 2018) | 0.62
Self Training (Yun et al. 2020) | 0.61
Messy Collab (Arani, Sarfraz, and Zonooz 2019) | 0.60
Noisy Teacher (Sau and Balasubramanian 2016) | 0.59
TAKD (Mirzadeh et al. 2019) | 0.59
RCO (Jin et al. 2019) | 0.58
Probability Shift (Wen, Lai, and Qian 2019) | 0.58
Table 2:  The accuracies of networks trained by some of various knowledge distillation algorithms KD-Lib packages on the CIFAR10 dataset. All models were trained with the same hyperparameter set to ensure a fair comparison. We consider ResNet34 as the teacher network (with an accuracy of 0.63) and report accuracies for the student network (ResNet18). _None_ refers to a ResNet18 model trained from scratch without any model compression algorithm. The compression ratio for all of the knowledge distillation algorithms is 50.7%  Pruning Epoch | % Model Pruned | Accuracy
---|---|---
1 | 0.0 | 0.9878
2 | 0.10 | 0.9891
3 | 0.19 | 0.9890
Table 3: Pruning percentage and accuracy of ResNet18 model on MNIST using Lottery Ticket Pruning (Frankle and Carbin 2018). Each pruning epoch consists of 5 training epochs. ’Model pruned’ is the percentage of model pruned and ’Accuracy’ is the corresponding accuracy at the end of the epoch. Algorithm | % Size Change | BA | NA
---|---|---|---
Static | -0.75 | 0.72 | 0.70
QAT | -0.75 | 0.72 | 0.71
Dynamic | -0.19 | 0.70 | 0.70
Table 4: Comparison of various quantization algorithms. ’BA’ (Base Accuracy)
is the accuracy of the model before quantization, and ’NA’ (New Accuracy) is
the accuracy of the model after quantization. ’% Size change’ refers to the
change in size after quantization. In Static Quantization and QAT, ResNet18 is
tested on the CIFAR10 dataset. For Dynamic Quantization, LSTM is tested on
IMDB dataset.

## Conclusion and Future Work

In this paper, we present KD-Lib, an easy-to-use PyTorch-based library for
Knowledge Distillation, Pruning and Quantization. KD-Lib is designed to
facilitate the adoption of current model compression techniques and act as a
catalyst for further research in this direction. We plan on actively
maintaining the library and also expanding it to include more algorithms and
desirable features (distributed training, for example) in the future. We
further plan on extending this library to other domains relevant to the
research community including but not limited to explainability and
interpretability in knowledge distillation.

> ## References
>
>   * Arani, Sarfraz, and Zonooz (2019) Arani, E.; Sarfraz, F.; and Zonooz, B.
> 2019.  Improving Generalization and Robustness with Noisy Collaboration in
> Knowledge Distillation.
>   * Frankle and Carbin (2018) Frankle, J.; and Carbin, M. 2018.  The Lottery
> Ticket Hypothesis: Finding Sparse, Trainable Neural Networks.
>   * Hinton, Vinyals, and Dean (2015) Hinton, G.; Vinyals, O.; and Dean, J.
> 2015.  Distilling the Knowledge in a Neural Network.
>   * Jacob et al. (2018) Jacob, B.; Kligys, S.; Chen, B.; Zhu, M.; Tang, M.;
> Howard, A.; Adam, H.; and Kalenichenko, D. 2018.  Quantization and Training
> of Neural Networks for Efficient Integer-Arithmetic-Only Inference.  In
> _2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition_.
>   * Jin et al. (2019) Jin, X.; Peng, B.; Wu, Y.; Liu, Y.; Liu, J.; Liang,
> D.; Yan, J.; and Hu, X. 2019\.  Knowledge Distillation via Route Constrained
> Optimization.
>   * LeCun et al. (1990) LeCun, Y.; Denker, J. S.; ; and Solla, S. A. 1990.
> Optimal brain damage.  _In Advances in Neural Information Processing
> Systems_ .
>   * Mirzadeh et al. (2019) Mirzadeh, S.-I.; Farajtabar, M.; Li, A.; Levine,
> N.; Matsukawa, A.; and Ghasemzadeh, H. 2019.  Improved Knowledge
> Distillation via Teacher Assistant.
>   * Sau and Balasubramanian (2016) Sau, B. B.; and Balasubramanian, V. N.
> 2016.  Deep Model Compression: Distilling Knowledge from Noisy Teachers.
>   * Tang et al. (2019) Tang, R.; Lu, Y.; Liu, L.; Mou, L.; Vechtomova, O.;
> and Lin, J. 2019.  Distilling Task-Specific Knowledge from BERT into Simple
> Neural Networks.  _CoRR_ abs/1903.12136.
>   * Wen, Lai, and Qian (2019) Wen, T.; Lai, S.; and Qian, X. 2019.
> Preparing Lessons: Improve Knowledge Distillation with Better Supervision.
>   * Yun et al. (2020) Yun, S.; Park, J.; Lee, K.; and Shin, J. 2020.
> Regularizing Class-Wise Predictions via Self-Knowledge Distillation.  In
> _2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition_.
>   * Zhang et al. (2018) Zhang, Y.; Xiang, T.; Hospedales, T. M.; and Lu, H.
> 2018.  Deep Mutual Learning.  In _2018 IEEE/CVF Conference on Computer
> Vision and Pattern Recognition_.
>   * Zmora et al. (2019) Zmora, N.; Jacob, G.; Zlotnik, L.; Elharar, B.; and
> Novik, G. 2019.  Neural Network Distiller: A Python Package For DNN
> Compression Research URL https://arxiv.org/abs/1910.12232.
>