# Distilling Knowledge from Reader to Retriever for Question Answering

Gautier Izacard1,2,3, Edouard Grave1
1Facebook AI Research, 2École normale supérieure, PSL University, 3Inria
<EMAIL_ADDRESS>
###### Abstract

The task of information retrieval is an important component of many natural
language processing systems, such as open domain question answering. While
traditional methods were based on hand-crafted features, continuous
representations based on neural networks recently obtained competitive
results. A challenge of using such methods is to obtain supervised data to
train the retriever model, corresponding to pairs of query and support
documents. In this paper, we propose a technique to learn retriever models for
downstream tasks, inspired by knowledge distillation, and which does not
require annotated pairs of query and documents. Our approach leverages
attention scores of a reader model, used to solve the task based on retrieved
documents, to obtain synthetic labels for the retriever. We evaluate our
method on question answering, obtaining state-of-the-art results.

##  1 Introduction

Information retrieval is an important component for many natural language
processing tasks, such as question answering (Voorhees et al., 1999) or fact
checking (Thorne et al., 2018). For example, many real world question
answering systems start by retrieving a set of support documents from a large
source of knowledge such as Wikipedia. Then, a finer-grained model processes
these documents to extract the answer. Traditionally, information retrieval
systems were based on hand-crafted sparse representations of text documents,
such as TF-IDF or BM25 (Jones, 1972; Robertson et al., 1995). Recently,
methods based on dense vectors and machine learning have shown promising
results (Karpukhin et al., 2020; Khattab et al., 2020). Deep neural networks
based on pre-training, such as BERT (Devlin et al., 2019), have been used to
encode documents into fixed-size representations. These representations are
then queried using approximate nearest neighbors (Johnson et al., 2019). These
techniques have lead to improved performance on various question answering
tasks.

A challenge of applying machine learning to information retrieval is to obtain
training data for the retriever. To train such models, one needs pairs of
queries and the corresponding list of documents that contains the information
corresponding to the queries. Unfortunately, hand-labeling data to that end is
time consuming, and many datasets and applications lack such annotations. An
alternative approach is to resort to heuristics, or weakly supervised
learning, for example by considering that all documents containing the answer
are positive examples. However, these approaches suffer from the following
limitations. First, frequent answers or entities might lead to false positive
examples. As an example, consider the question _“where was Ada Lovelace
born?”_. The sentence _“Ada Lovelace died in 1852 in London”_ would be
considered as a positive example, because it contains the answer _“London”_. A
second limitation is that for some tasks, such as fact checking or long form
question answering, such heuristics might not be applicable directly.

In this paper, we propose a procedure to learn retriever systems without
strong supervision in the form of pairs of queries and documents. Following
previous work (Chen et al., 2017), our approach uses two models: the first one
retrieves documents from a large source of knowledge (_the retriever_), the
second one processes the support documents to solve the task (_the reader_).
Our method is inspired by _knowledge distillation_ , and uses the reader model
to obtain synthetic labels to train the retriever model. More precisely, we
use a sequence-to-sequence model as the reader, and use the attention
activations over the input documents as synthetic labels to train the
retriever. Said otherwise, we assume that attention activations are a good
proxy for the relevance of documents. We then train the retriever to reproduce
the ranking of documents corresponding to that metric.

We make the following contributions:

  * •

First, we show that attention scores from a sequence-to-sequence reader model
are a good measure of document relevance (Sec. 3.2) ;

  * •

Second, inspired by knowledge distillation, we propose to iteratively train
the retriever from these activations, and compare different loss functions
(Sec. 3.4) ;

  * •

Finally, we evaluate our method on three question-answering benchmarks,
obtaining state-of-the-art results (Sec. 4).

##  2 Related work

We briefly review information retrieval based on machine learning. We refer
the reader to Manning et al. (2008) and Mitra et al. (2018) for a more
exhaustive introduction to the subject.

#### Vector space models.

In traditional information retrieval systems, documents and queries are
represented as sparse vectors, each dimension corresponding to a different
term. Different schemes have been considered to weigh the different term, the
most well known being based on inverse document frequency, or _term
specificity_ (Jones, 1972). This technique was later extended, leading to the
BM25 weighting scheme which is still widely used today (Robertson et al.,
1995). A limitation of sparse representations is that the terms of the query
need to match the terms of the returned documents. To overcome this,
Deerwester et al. (1990) proposed to use latent semantic analysis for
indexing, leading to low-dimension dense representations of documents.

#### Neural information retrieval.

Following the success of deep learning for other natural processing tasks,
neural networks were applied to the task of information retrieval. Huang et
al. (2013) proposed a deep bag-of-words model, where queries and documents
were embedded independently, a technique known as _bi-encoder._ Documents were
then ranked by using the cosine similarity with the query, and the model was
trained on clickthrough data from a search engine. This technique was later
extended by using convolutional neural networks (Shen et al., 2014) and
recurrent neural networks (Palangi et al., 2016). A limitation of
independently embedding documents and query is that it does not capture fine-
grained interactions between the query and documents. This lead Nogueira & Cho
(2019) and Yang et al. (2019) to use a BERT model to jointly embed documents
and query, a technique known as _cross-encoder._

#### End-to-end retrieval.

Most of the methods described in the previous paragraph were used to re-rank a
small number of documents, usually returned by a traditional IR systems. In
the context of _ad-hoc_ document retrieval, Gillick et al. (2018) showed that
bi-encoder models could be competitive with traditional IR systems. For open
domain question-answering, Karpukhin et al. (2020) introduced _dense passage
retrieval_ (DPR), which uses dense embeddings and nearest neighbors search.
More precisely, question and passage embeddings are obtained using a BERT-
based bi-encoder model, which is trained on a small dataset of question and
passage pairs. Then, the full knowledge source (Wikipedia) is encoded using
this model, and passages are queried by computing the k-nearest neighbors of
the embedding of the question. Jointly embedding the query and documents makes
the application of cross-encoder models intractable to large database. To
address this limitation, Humeau et al. (2019) introduced the _poly-encoder_
architecture, in which each documents is represented by multiple vectors
instead of one. Similarly, Khattab et al. (2020) proposed a scoring function
where each term of the query and documents is represented by a single vector.
To make the method tractable, their system retrieves documents with an
approximate score, which are then re-ranked with the exact one. Finally, Luan
et al. (2020) conducts a theoretical and empirical study of sparse, dense and
cross-attention information retrieval systems.

#### Unsupervised learning.

Closest to our work, there is growing body of work trying to learn information
retrieval systems from unsupervised data. Lee et al. (2019) introduced the
_inverse cloze task_ for pre-training retrievers, which can then be fine-tuned
end-to-end on question-answering tasks. This pre-training scheme was later
evaluated for ad-hoc document retrieval by Chang et al. (2020). Guu et al.
(2020) proposed to augment language model pre-training with a retriever
module, which is trained using the masked language modeling objective.
Similarly, Lewis et al. (2020a) introduced a sequence-to-sequence model that
is pre-trained by generating a target text, after retrieving a set of related
texts. Lewis et al. (2020b) further train the retriever obtained in Karpukhin
et al. (2020) by backpropagating to the retriever the error between the
generated output and the gold answer.

Simultaneously to our work, Yang & Seo (2020) proposes to train a retriever
with knowledge distillation. The main difference with our method is the nature
of the synthetic labels that are used to train the retriever. Yang & Seo
(2020) uses the DPR reader, which includes a classifier that predicts which
passage contains the answer, and can be seen as a cross-encoder reranker. This
technique thus consists in distilling a cross-encoder retriever to a bi-
encoder retriever. In contrast, our method uses the internal attention scores
of the reader, which does not require additional supervision besides pairs of
question and answer.

##  3 Methodology

Our system is composed of two modules, the retriever and the reader, following
the standard pipeline for open-domain question answering. Given an input
question these modules are used in a two-step process to generate an answer.
First the retriever selects support passages in a large knowledge source. Then
these passages are processed by the reader, along with the question, to
generate an answer. For the reader module we use the Fusion-in-Decoder model
(Izacard & Grave, 2020), which achieves state-of-the-art performance when
combined with BM25 or DPR (Karpukhin et al., 2020). It is based on a sequence-
to-sequence architecture, and is initialized from pre-trained models such as
T5 or BART (Raffel et al., 2019; Lewis et al., 2019).

The focus of this work is to train the retriever without strong supervision or
weakly supervised learning based on heuristics. For this we propose to train
the retriever by learning to approximate the attention score of the reader.
The training scheme outlined here can be seen as a student-teacher pipeline,
where the teacher, the reader module, produces targets which are used to train
a student network, the reader. By doing so, we hope to leverage the signal
extracted from the question-answer pairs by the reader. Since the goal of the
retriever is to retrieve the most relevant passages, by training the retriever
to estimate the reader attention scores, we implicitly make the assumption
that these scores are a good proxy for the usefulness of a passage to answer
the question.

In this section we will first describe the Fusion-in-Decoder architecture,
before elaborating on the signal which is used to train the retriever, the
design of the retriever, and how this module is trained.

###  3.1 Cross-Attention Mechanism

First, let us briefly review the Fusion-in-Decoder model (FiD, Izacard &
Grave, 2020). The underlying architecture is a sequence-to-sequence model,
composed of an encoder and a decoder. The encoder independently processes
$n_{p}$ different text inputs $(s_{k})_{1\leq k\leq n_{p}}$. In the case of
open-domain question answering based on Wikipedia, each input $s_{k}$ is the
concatenation of the question $q$ and a support passage, with special tokens
question:, title: and context: added before the question, the title of the
Wikipedia article and the text of each passage. The output representations of
the encoder are then concatenated to form a global representation $\mathbf{X}$
of dimension $(\sum_{k}\ell_{k})\times d$, where $\ell_{k}$ is the length of
the $k$-th segment and $d$ is the dimension of the embeddings and hidden
representations of the model. Then, the decoder processes this representation
as a regular autoregressive model, alternating self-attention, cross-attention
and feed-forward modules.

Only the cross-attention module explicitly takes as input the global output
representation $\mathbf{X}$ of the encoder. If $\mathbf{H}\in\mathbb{R}^{d}$
denotes the output of the previous self-attention layer of the decoder, the
cross-attention operation consists in the following operations. First, queries
$\mathbf{Q}$, keys $\mathbf{K}$ and values $\mathbf{V}$ are computed by
applying linear transformations:

$\displaystyle\mathbf{Q}=\mathbf{W}_{Q}\mathbf{H},\quad\mathbf{K}=\mathbf{W}_{K}\mathbf{X},\quad\mathbf{V}=\mathbf{W}_{V}\mathbf{X}.$

Then a similarity score between the query at position $i$, $\mathbf{Q}_{i}$,
and the key at position $j$, $\mathbf{K}_{j}$, is obtained by computing the
dot-product between these two elements, and normalized over the dimension:

$\displaystyle\alpha_{i,j}=\mathbf{Q}_{i}^{T}\mathbf{K}_{j},\qquad\tilde{\alpha}_{i,j}=\frac{\exp(\alpha_{i,j})}{\sum_{m}\exp(\alpha_{i,m})}.$

A new representation is obtained as a sum of the values, weighted by the
attention probabilities, before going through a final linear transformation
$\mathbf{W}_{o}$:

$\displaystyle\mathbf{O}_{i}=\mathbf{W}_{O}\sum_{j}\tilde{\alpha}_{i,j}\mathbf{V}_{i,j}$

The operations described above are performed in parallel with different linear
transformations in the case of multi-head attention. Finally a normalization
layer is applied, and this pipeline is wrapped by a skip connection. See
Vaswani et al. (2017) for more details on the structure of Transformers.

###  3.2 Cross-Attention Score as a Relevance Measure for Passage Retrieval

In some sense, the attention scores $\alpha_{:,j}$ involving the $j$-th key
measures the importance of this key, and corresponding value, to compute the
next representation. We hypothesize that it is good proxy to estimate the
relevance of a passage — the more the tokens in a text segment are attended
to, the more relevant the text segment is to answer the question.

Given the reader model, an input question $q$ and a corresponding set of
support passages ${\mathcal{D}}_{q}~{}=~{}(p_{k})_{1\leq k\leq n}$, we obtain
relevance scores $(G_{q,p_{k}})_{1\leq k\leq n}$ for each passage by
aggregating attention scores. In particular, the score $G_{q,p_{k}}$ is
obtained by averaging the pre-attention scores $\alpha_{0,:}$ over all the
tokens in the input $s_{k}$ corresponding to the passage $p_{k}$, all the
layers and all the heads of the decoder. Note that the FiD decoder jointly
processes the passages, and thus the score $G_{q,p_{k}}$ depends on the other
support passages. We consider other pooling operators, such as _max_ , to
aggregate attention scores over layers, heads and tokens and empirically
compare them in Sec. 5.2.

Before we proceed, let us consider the following simple experiment, which is a
first indication that reader attention scores are indeed a strong relevance
signal. Given a question and 100 passages retrieved with DPR, our goal is to
select the 10 best passages. When using the top 10 passages from DPR instead
of the top 100, the performance of our reader drops from 48.2 EM to 42.9 EM.
On the other hand, if we select the top 10 documents according to the
attention scores, the performance only drops to 46.8 EM.

###  3.3 Dense bi-encoder for passage retrieval

Ideally, we would like to rank passages according to the reader cross-
attention scores. In practice however, since the passages and the question
need to be processed simultaneously by the reader module it is impractical to
query a large knowledge source this way. Thus, we use a retriever model
composed of an embedder function $E$ that maps any text passage to a
$d$-dimensional vector, such that the similarity score between a question $q$
and a passage $p$ is defined as $S_{\theta}(q,p)=E(q)^{T}E(p)$. This
similarity metric enables us to index all passages in the knowledge source as
a preprocessing step. Then at runtime, passages with the highest similarity
score with the input question are retrieved, by using an efficient similarity
search library such as FAISS (Johnson et al., 2019).

For the embedder we use BERT and follow DPR by considering that the encodings
$E(q)$ and $E(p)$ are obtained by extracting the representation of the initial
[CLS] token. This leads to a representation of dimension $d=768$ in the case
of a base model. Differently from DPR, we use the same encoding function $E$
for the questions and passages by sharing parameters.

###  3.4 Distilling the cross-attention score to a bi-encoder

In this section, we describe how to train the retriever model, based on the
relevance scores obtained in Sec. 3.2. For the training objective of the
retriever, we propose to minimize the KL-divergence between the output
$S_{\theta}(q,p)$ and the score $G_{q,p}$ after normalization:

$\displaystyle{\mathcal{L}}_{\text{KL}}(\theta,{\mathcal{Q}})=\sum_{q\in{\mathcal{Q}},p\in{\mathcal{D}}_{q}}\tilde{G}_{q,p}(\log\tilde{G}_{q,p}-\log\tilde{S}_{\theta}(q,p)),$

where

$\displaystyle\tilde{G}_{q,p}=\frac{\exp(G_{q,p})}{\sum_{p^{\prime}\in{\mathcal{D}}_{q}}\exp(G_{q,p^{\prime}})},\qquad\tilde{S}_{\theta}(q,p)=\frac{\exp(S_{\theta}(q,p))}{\sum_{p^{\prime}\in{\mathcal{D}}_{q}}\exp(S_{\theta}(q,p^{\prime}))}.$

In Sec. 5.1 we present results obtained when using alternatives to this
training objective. We consider two other objectives which have been used in
Dehghani et al. (2017), where BM25 is used as a teacher model to train a
neural ranker. A first option consists in training the retriever with a
regression approach by minimizing the mean squared error:

$\displaystyle{\mathcal{L}}_{\text{MSE}}(\theta,{\mathcal{Q}})=\sum_{q\in{\mathcal{Q}},p\in{\mathcal{D}}_{q}}(S_{\theta}(q,p)-G_{q,p})^{2}.$

The second option we consider is to use a max-margin loss that explicitly
penalizes inversions in the ranking estimated by the retriever:

$\displaystyle{\mathcal{L}}_{\text{ranking}}(\theta,{\mathcal{Q}})=\sum_{q\in{\mathcal{Q}},p_{1},p_{2}\in{\mathcal{D}}_{q}}\max\left(0,\gamma-\text{sign}(G_{q,p_{1}}-G_{q,p_{2}})(S_{\theta}(q,p_{1})-S_{\theta}(q,p_{2}))\right).$

In words, if $p_{1}$ is more relevant to answer the question $q$ than $p_{2}$
, i.e. $G_{q,p_{1}}>G_{q,p_{2}}$, the loss pushes the retriever score of
$p_{1}$ to be larger than the score of $p_{2}$ by at least a margin of
$\gamma$.

###  3.5 Iterative training

In this section, we explain how iterative training can be used with the
student-teacher scheme described in the previous section, similarly to Khattab
et al. (2020). This iterative procedure can be interpreted as using the
current retriever to sample negative examples, in order to train a new
retriever. When learning a retriever with discriminative training, negative
samples play an important role, and various strategies have been considered in
previous work. Karpukhin et al. (2020) compared random sampling with using the
top-k passages from BM25 which do not contain the answer and with using the
positive passages from other queries. Consider that for each question, we have
an initial set of support documents ${\mathcal{D}}_{q}^{0}$. We propose to use
an iterative pipeline where each iteration can be described as the following
4-step process:

  1. 1.

Train the reader $R$ using the set of support documents for each question
${\mathcal{D}}_{q}^{0}$.

  2. 2.

Compute aggregated attention scores
$(G_{q,p})_{q\in{\mathcal{Q}},p\in{\mathcal{D}}^{0}_{q}}$ with the reader $R$.

  3. 3.

Train the retriever $E$ using the scores
$(G_{q,p})_{q\in{\mathcal{Q}},p\in{\mathcal{D}}^{0}_{q}}$.

  4. 4.

Retrieve top-passages with the new trained retriever $E$.

This multi-step procedure can be repeated multiple times. A critical point of
the training procedure is the initial set of documents corresponding to each
question. In Sec. 4, we compare retrievers obtained by starting from documents
obtained using BM25 or cosine similarity from a BERT model. In particular, we
show that while the initial performance with BERT is low, the iterative
procedure allows to greatly improve the performance of the model.

##  4 Experiments

In this section we evaluate the student-teacher training procedure from the
previous section. We show that we obtain competitive performance without
strong supervision for support documents.

###  4.1 Experimental Setting

#### Datasets.

We perform experiments on TriviaQA (Joshi et al., 2017) and NaturalQuestions
(Kwiatkowski et al., 2019), two standard benchmarks for open-domain question
answering. TriviaQA is made of questions from trivia and quiz league websites,
and does not contain gold support documents. NaturalQuestions contains
questions corresponding to web search queries, and gold support documents from
Wikipedia. Following the setting from Lee et al. (2019); Karpukhin et al.
(2020), we use the original evaluation set as test set, and keep 10% of the
training data for validation. We use the Wikipedia dump from Dec. 20, 2018 for
support documents, splitting articles into non-overlapping passages of 100
tokens, and applying the same preprocessing as Chen et al. (2017).

We also evaluate on NarrativeQuestions (Kočiskỳ et al., 2018), using a
publicly available preprocessed
version.111https://cs.nyu.edu/~kcho/NarrativeQA This is a reading
comprehension dataset built on a corpus of books and movie scripts. For each
story, questions are generated by human annotators based on a summary of the
given document. We consider the full story setting, where the task is to
answer questions given the entire story and not the summary used to generate
question-answer pairs. Here the knowledge source is not the same for all
questions: given a question the retrieval operation is performed on all
passages of the associated story. These passages are obtained by dividing the
story in chunks of 100 words. These stories are long documents, with an
average of 60k words. While part of the documents could be processed entirely
by the Fusion-in-Decoder module, it is interesting to limit the number of
support passages to reduce the computational cost of the reading step.

While answers in TriviaQA and NaturalQuestions are short, NarrativeQA answers
are about five words long on average, with medium length answers such as _”He
dismantles it and attaches it to his mother’s jeep”_ which answers the
question _”What does Mark do with his radio station?”_. Notably a significant
number of answers do not correspond to spans in the story. It is thus not
straightforward to train the retriever with heuristics using question-answer
pairs. In our case we use the same pipeline as for TriviaQA and
NaturalQuestions, demonstrating the flexibility of our approach.

#### Evaluation.

The model performance is assessed in two ways. First, following previous work
such as DPR and ColbertQA, we report the top-$k$ retrieval accuracy (P@k),
which is the percentage of questions for which at least one passage of the
top-$k$ retrieved passages contains the gold answer. It is unclear how well
this metric evaluates the retriever performance, since the answer can be
contained in a passage without being related to the question. This is notably
true for common words or entities.

We also report the final end-to-end performance of the question answering
system composed of the retriever and reader modules. This is the metric we are
fundamentally interested in. For TriviaQA and NaturalQuestions, predicted
answers are evaluated with the standard exact match metric (EM), as introduced
by Rajpurkar et al. (2016). For NarrativeQA we report the metrics proposed in
the original paper: ROUGE-L, BLEU-1, BLEU-4 and METEOR.

|  |  | BERT |  | BM25
---|---|---|---|---|---
| Iter. |  | P@20 | P@100 | Dev EM |  | P@20 | P@100 | Dev EM
NaturalQuestions | 0 |  | 4.8 | 12.0 | 9.8 |  | 59.3 | 74.0 | 41.2
1 |  | 32.2 | 45.8 | 16.9 |  | 76.4 | 84.3 | 46.8
2 |  | 51.1 | 62.6 | 28.6 |  | 80.4 | 86.7 | 47.9
3 |  | 67.8 | 76.8 | 39.3 |  | 80.0 | 86.3 | 46.2
TriviaQA | 0 |  | 4.6 | 12.0 | 9.7 |  | 75.0 | 82.3 | 65.3
1 |  | 37.1 | 59.4 | 19.6 |  | 79.0 | 85.5 | 66.7
2 |  | 60.8 | 73.4 | 43.3 |  | 82.1 | 86.5 | 67.5
3 |  | 72.0 | 83.2 | 52.0 |  | 81.6 | 86.6 | 67.7
4 |  | 76.4 | 84.6 | 62.3 |  | - | - | -
Table 1: Iterative training starting with documents retrieved with BERT and
BM25. Iteration 0 corresponds to the performance of the reader trained on the
set of initial support documents. We report all metrics on the validation set.

###  4.2 Technical details

#### Initialization.

Similarly to DPR, we initialize the retriever with the BERT base model,
pretrained with uncased text. The Fusion-in-Decoder reader is initialized with
the T5 base model. A critical component of the iterative training procedure is
the initialization of the support passages ${\mathcal{D}}_{q}^{0}$ associated
with each question $q$. For this we consider different options. The first one
is to use passages retrieved using BM25. We use the implementation from Apache
Lucene222lucene.apache.org with default parameters, and tokenize questions and
passages with SpaCy333spacy.io. We also use passages obtained with BERT as a
retriever without fine-tuning, this leads to poor initial performance. Finally
in Table 2 we show that initializing ${\mathcal{D}}_{q}^{0}$ with passages
obtained with DPR (Karpukhin et al., 2020) outperforms the two previous
initializations. We train all retrievers using 100 passages. For the reader,
we use 100 passages for NaturalQuestions and TriviaQA and 20 passages for
NarrativeQA.

#### Iterative training.

We apply the iterative training procedure on each dataset independently. Both
the reader and the retriever are fine-tuned using the ADAM algorithm (Kingma &
Ba, 2014), with a batch of size 64. The reader is trained for 10k gradient
steps with a constant learning rate of $10^{-4}$, and the best model is
selected based on the validation performance. The retriever is trained with a
constant learning rate of $5\cdot 10^{-5}$ until the performance saturates. To
monitor the performance of the retriever during training, we measure the
similarity between the reader and the retriever rankings. At each new training
iteration the reader is reinitialized from T5 base, while we pursue the
training of the retriever. We found that restarting from T5 base is important
for the first iterations when starting with BERT documents. We have not tried
to reinitialize the retriever between each iteration. More details on the
hyperparameters and the training procedure are reported in Appendix A.2.

###  4.3 Results

In Table 1, we report the performance of our approach for different number of
self-training iterations. Generally, we observe that the accuracy of our
system increases with the number of iterations, obtaining strong performance
after a few iterations. Interestingly, while the initial performance with
documents retrieved with BERT is very poor, our method still reach competitive
scores on TriviaQA, and to a lesser extent, NaturalQuestions. However, a
second observation is that the quality of the initial document sets plays an
important role on the performance of the end system. Indeed, we observe that
starting the procedure from BM25 documents, which are higher quality as
indicated by the performance of the system at iteration 0, leads to stronger
results than using BERT documents. An interesting research question would be
to explore pre-training of the initial BERT model for retrieval, for example
by using the inverse cloze task.

In Table 2, we report the performance of our approach, as well as existing
state-of-the-art systems on TriviaQA and NaturalQuestions. In addition to
initializing our method with documents retrieved with BM25 and BERT, we also
train a system by starting from DPR documents. First, we observe that our
method improve the performance over the state-of-the-art, even when starting
from BM25 documents. This validates our assumption that it is possible to
obtain strong retrievers without the need of supervision for the documents.
Second, when starting from DPR passages, our method leads to a +3.9 EM
improvement on TriviaQA and +1.1 EM improvement on NaturalQuestions when the
final evaluation is carried out with a large reader.

In Table 3, we report the performance of our method on the NarrativeQA
dataset. We use the setting where the knowledge source corresponds to the
whole document, and in particular, we do not use the summary. We compare our
results to the best ones reported in the original paper for this setting.
Similar to results obtained on NaturalQuestions and TriviaQA, we observe that
training the retriever by using the attention scores of the reader leads to
improvements, compared to the BM25 baseline.

Model | NQ | TriviaQA
---|---|---
| dev. | test | dev. | test
DPR (Karpukhin et al., 2020) | - | 41.5 | - | 57.9
RAG (Lewis et al., 2020b) | - | 44.5 | - | 56.1
ColBERT-QA (Khattab et al., 2020) | - | 48.2 | - | 63.2
Fusion-in-Decoder (T5 base) (Izacard & Grave, 2020) | - | 48.2 | - | 65.0
Fusion-in-Decoder (T5 large) (Izacard & Grave, 2020) | - | 51.4 | - | 67.6
Ours (starting from BERT, T5 base) | 39.3 | 40.0 | 62.5 | 62.7
Ours (starting from BM25, T5 base) | 47.9 | 48.9 | 67.7 | 67.7
Ours (starting from DPR, T5 base) | 48.0 | 49.6 | 68.6 | 68.8
Ours (starting from DPR, T5 large) | 51.9 | 53.7 | 71.9 | 72.1
Table 2: Comparison to state-of-the-art models on NaturalQuestions and TriviaQA. Method | Iter. | Rouge-L | Bleu-1 | Bleu-4 | Meteor
---|---|---|---|---|---
|  | dev. | test | dev. | test | dev. | test | dev. | test
Best from Kočiskỳ et al. (2018) | - | 14.5 | 14.0 | 20.0 | 19.1 | 2.23 | 2.1 | 4.6 | 4.4
DPR + FiD | - | 29.7 | 30.8 | 33.0 | 34.0 | 6.7 | 6.9 | 10.3 | 10.8
Ours starting from BM25 | 0 | 29.9 | 30.3 | 34.6 | 33.7 | 7.1 | 6.5 | 10.5 | 10.4
Ours starting from BM25 | 1 | 31.6 | 32.0 | 34.9 | 35.3 | 7.6 | 7.5 | 11.0 | 11.1
Table 3: Performance on NarrativeQA.

##  5 Ablations

In this section, we investigate design choices regarding two key elements of
our approach: the training objective and the aggregation of cross-attention
scores. For all experiments, we consider a simplified experimental setting: a
single training iteration is performed on NaturalQuestions, starting from BM25
passages.

###  5.1 Training objectives

In Table 4 we report the performance of our model trained with the different
training objectives described in Sec. 3.3. We observe that using the KL-
divergence between the aggregated scores of the reader and the scores of the
retriever outperforms the other objective functions.

Method | P@5 | P@20 | P@100 | Dev EM
---|---|---|---|---
Mean Squared Error | 46.5 | 61.2 | 73.9 | 40.6
Max-margin loss, $\gamma=1$ | 60.3 | 73.6 | 82.7 | 45.4
Max-margin loss, $\gamma=0.2$ | 60.3 | 73.5 | 82.6 | 45.8
Max-margin loss, $\gamma=0.1$ | 60.2 | 73.5 | 82.6 | 45.1
KL-divergence | 64.7 | 76.4 | 84.3 | 46.8
Table 4: Comparison of training objectives on NaturalQuestions after one
iteration. We report all the metrics on the validation set.

###  5.2 How to aggregate cross-attention scores?

In Section 4 the cross-attention scores $\alpha$ are aggregated in a specific
way, in order to obtain a single scalar used to train the retriever. Formally
let us denote by $\alpha_{i,j,k,h}$ the cross-attention scores between token
$i$ of the output and token $j$ of the input, for the $k$-th layer and $h$-th
head. Then, the scores $G_{q,p}$ for $p\in{\mathcal{D}}_{q}$ used in Section 4
are computed as follows:

$G_{q,p}=\operatorname*{mean}_{j,k,h}\alpha_{0,j,k,h},$

where $j$ describes the input tokens corresponding to $p$. In Table 5 we
explore alternatives to this choice by considering different aggregation
schemes. In particular, we consider (1) taking the max over the input tokens
corresponding to passage $p$ instead of the average, (2) taking the average
over the output tokens instead of taking the score of the first token, (3)
taking the mean over the last six layers instead of all the layers, (4) taking
the max over the layers instead of the average, (5) taking the max over the
heads instead of the average. We observe that the performance of our approach
is relatively stable to the choice of aggregation, and that the best result is
obtained by averaging, except over the output tokens where it is best to only
consider the first token.

Method | P@5 | P@20 | P@100 | Dev EM
---|---|---|---|---
(0) $\operatorname*{mean}_{j,k,h}\alpha_{0,j,k,h}$ | 64.7 | 76.4 | 84.3 | 46.8
(1) $\operatorname*{mean}_{k,h}\max_{j}\alpha_{0,j,k,h}$ | 61.2 | 72.5 | 81.0 | 46.0
(2) $\operatorname*{mean}_{i,j,k,h}\alpha_{i,j,k,h}$ | 63.5 | 75.3 | 83.1 | 45.8
(3) $\operatorname*{mean}_{7\leq k\leq 12,j,h}\alpha_{0,j,k,h}$ | 64.1 | 75.7 | 83.8 | 46.4
(4) $\operatorname*{mean}_{j,h}\max_{k}\alpha_{0,j,k,h}$ | 63.9 | 75.5 | 83.7 | 46.5
(5) $\operatorname*{mean}_{j,k}\max_{h}\alpha_{0,j,k,h}$ | 64.2 | 76.1 | 83.9 | 46.8
Table 5: Comparison of attention aggregation schemes on NaturalQuestions after
one iteration. The index $i$ corresponds to output tokens, $j$ corresponds to
input tokens of a given passage, $h$ to heads and $k$ to layers of the
decoder. We report all metrics on the validation set.

##  6 Conclusion

In this paper, we introduce a method to train an information retrieval module
for downstream tasks, without using pairs of queries and documents as
annotations. Our approach is inspired by knowledge distillation, where the
retriever module corresponds to the student model and the reader module
corresponds to the teacher model. In particular, we use the cross-attention
scores, from a sequence-to-sequence reader, to obtain synthetic targets for
the retriever. We compare different ways to aggregate the scores, as well as
different training objectives to learn the retriever. We show that iteratively
training the reader and the retriever leads to better performance, and obtain
state-of-the-art performance on competitive question answering benchmarks. In
the future, we would like to explore better pre-training strategies for the
retriever module, as well as better scoring functions for the retriever.

## References

  * Chang et al. (2020) Wei-Cheng Chang, Felix X Yu, Yin-Wen Chang, Yiming Yang, and Sanjiv Kumar.  Pre-training tasks for embedding-based large-scale retrieval.  _arXiv preprint arXiv:2002.03932_ , 2020.
  * Chen et al. (2017) Danqi Chen, Adam Fisch, Jason Weston, and Antoine Bordes.  Reading Wikipedia to answer open-domain questions.  In _Proc. ACL_ , 2017.
  * Deerwester et al. (1990) Scott Deerwester, Susan T Dumais, George W Furnas, Thomas K Landauer, and Richard Harshman.  Indexing by latent semantic analysis.  _Journal of the American society for information science_ , 41(6):391–407, 1990.
  * Dehghani et al. (2017) Mostafa Dehghani, Hamed Zamani, Aliaksei Severyn, Jaap Kamps, and W. Bruce Croft.  Neural ranking models with weak supervision.  In _Proceedings of the 40th International ACM SIGIR Conference on Research and Development in Information Retrieval_ , SIGIR ’17, pp. 65–74, New York, NY, USA, 2017. Association for Computing Machinery.  doi: 10.1145/3077136.3080832.  URL https://doi.org/10.1145/3077136.3080832.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In _Proc. NAACL_ , 2019.
  * Gillick et al. (2018) Daniel Gillick, Alessandro Presta, and Gaurav Singh Tomar.  End-to-end retrieval in continuous space.  _arXiv preprint arXiv:1811.08008_ , 2018.
  * Guu et al. (2020) Kelvin Guu, Kenton Lee, Zora Tung, Panupong Pasupat, and Ming-Wei Chang.  Realm: Retrieval-augmented language model pre-training.  _arXiv preprint arXiv:2002.08909_ , 2020.
  * Huang et al. (2013) Po-Sen Huang, Xiaodong He, Jianfeng Gao, Li Deng, Alex Acero, and Larry Heck.  Learning deep structured semantic models for web search using clickthrough data.  In _Proceedings of the 22nd ACM international conference on Information & Knowledge Management_, pp. 2333–2338, 2013.
  * Humeau et al. (2019) Samuel Humeau, Kurt Shuster, Marie-Anne Lachaux, and Jason Weston.  Poly-encoders: Transformer architectures and pre-training strategies for fast and accurate multi-sentence scoring.  _arXiv preprint arXiv:1905.01969_ , 2019.
  * Izacard & Grave (2020) Gautier Izacard and Edouard Grave.  Leveraging passage retrieval with generative models for open domain question answering.  _arXiv preprint arXiv:2007.01282_ , 2020.
  * Johnson et al. (2019) Jeff Johnson, Matthijs Douze, and Hervé Jégou.  Billion-scale similarity search with gpus.  _IEEE Transactions on Big Data_ , 2019.
  * Jones (1972) Karen Sparck Jones.  A statistical interpretation of term specificity and its application in retrieval.  _Journal of documentation_ , 1972.
  * Joshi et al. (2017) Mandar Joshi, Eunsol Choi, Daniel S. Weld, and Luke Zettlemoyer.  Triviaqa: A large scale distantly supervised challenge dataset for reading comprehension.  In _Proc. ACL_ , 2017.
  * Karpukhin et al. (2020) Vladimir Karpukhin, Barlas Oğuz, Sewon Min, Ledell Wu, Sergey Edunov, Danqi Chen, and Wen-tau Yih.  Dense passage retrieval for open-domain question answering.  _arXiv preprint arXiv:2004.04906_ , 2020.
  * Khattab et al. (2020) Omar Khattab, Christopher Potts, and Matei Zaharia.  Relevance-guided supervision for openqa with colbert, 2020.
  * Kingma & Ba (2014) Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Kočiskỳ et al. (2018) Tomáš Kočiskỳ, Jonathan Schwarz, Phil Blunsom, Chris Dyer, Karl Moritz Hermann, Gábor Melis, and Edward Grefenstette.  The NarrativeQA reading comprehension challenge.  _TACL_ , 2018.
  * Kwiatkowski et al. (2019) Tom Kwiatkowski, Jennimaria Palomaki, Olivia Redfield, Michael Collins, Ankur Parikh, Chris Alberti, Danielle Epstein, Illia Polosukhin, Matthew Kelcey, Jacob Devlin, Kenton Lee, Kristina N. Toutanova, Llion Jones, Ming-Wei Chang, Andrew Dai, Jakob Uszkoreit, Quoc Le, and Slav Petrov.  Natural Questions: a benchmark for question answering research.  _TACL_ , 2019.
  * Lee et al. (2019) Kenton Lee, Ming-Wei Chang, and Kristina Toutanova.  Latent retrieval for weakly supervised open domain question answering.  In _Proc. ACL_ , 2019.
  * Lewis et al. (2019) Mike Lewis, Yinhan Liu, Naman Goyal, Marjan Ghazvininejad, Abdelrahman Mohamed, Omer Levy, Ves Stoyanov, and Luke Zettlemoyer.  BART: Denoising sequence-to-sequence pre-training for natural language generation, translation, and comprehension.  _arXiv preprint arXiv:1910.13461_ , 2019.
  * Lewis et al. (2020a) Mike Lewis, Marjan Ghazvininejad, Gargi Ghosh, Armen Aghajanyan, Sida Wang, and Luke Zettlemoyer.  Pre-training via paraphrasing.  _arXiv preprint arXiv:2006.15020_ , 2020a.
  * Lewis et al. (2020b) Patrick Lewis, Ethan Perez, Aleksandara Piktus, Fabio Petroni, Vladimir Karpukhin, Naman Goyal, Heinrich Küttler, Mike Lewis, Wen-tau Yih, Tim Rocktäschel, et al.  Retrieval-augmented generation for knowledge-intensive nlp tasks.  _arXiv preprint arXiv:2005.11401_ , 2020b.
  * Loshchilov & Hutter (2019) Ilya Loshchilov and Frank Hutter.  Decoupled weight decay regularization.  In _International Conference on Learning Representations_ , 2019.
  * Luan et al. (2020) Yi Luan, Jacob Eisenstein, Kristina Toutanova, and Michael Collins.  Sparse, dense, and attentional representations for text retrieval.  _arXiv preprint arXiv:2005.00181_ , 2020.
  * Manning et al. (2008) Christopher D Manning, Hinrich Schütze, and Prabhakar Raghavan.  _Introduction to information retrieval_.  Cambridge university press, 2008.
  * Mitra et al. (2018) Bhaskar Mitra, Nick Craswell, et al.  An introduction to neural information retrieval.  _Foundations and Trends® in Information Retrieval_ , 13(1):1–126, 2018.
  * Nogueira & Cho (2019) Rodrigo Nogueira and Kyunghyun Cho.  Passage re-ranking with bert.  _arXiv preprint arXiv:1901.04085_ , 2019.
  * Palangi et al. (2016) Hamid Palangi, Li Deng, Yelong Shen, Jianfeng Gao, Xiaodong He, Jianshu Chen, Xinying Song, and Rabab Ward.  Deep sentence embedding using long short-term memory networks: Analysis and application to information retrieval.  _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , 24(4):694–707, 2016.
  * Raffel et al. (2019) Colin Raffel, Noam Shazeer, Adam Roberts, Katherine Lee, Sharan Narang, Michael Matena, Yanqi Zhou, Wei Li, and Peter J Liu.  Exploring the limits of transfer learning with a unified text-to-text transformer.  _arXiv preprint arXiv:1910.10683_ , 2019.
  * Rajpurkar et al. (2016) Pranav Rajpurkar, Jian Zhang, Konstantin Lopyrev, and Percy Liang.  SQuAD: 100,000+ questions for machine comprehension of text.  In _Proc. EMNLP_ , 2016.
  * Robertson et al. (1995) Stephen E Robertson, Steve Walker, Susan Jones, Micheline M Hancock-Beaulieu, Mike Gatford, et al.  Okapi at TREC-3.  _NIST Special Publication Sp_ , 1995.
  * Shen et al. (2014) Yelong Shen, Xiaodong He, Jianfeng Gao, Li Deng, and Grégoire Mesnil.  Learning semantic representations using convolutional neural networks for web search.  In _Proceedings of the 23rd international conference on world wide web_ , pp. 373–374, 2014.
  * Thorne et al. (2018) James Thorne, Andreas Vlachos, Christos Christodoulopoulos, and Arpit Mittal.  Fever: a large-scale dataset for fact extraction and verification.  _arXiv preprint arXiv:1803.05355_ , 2018.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Ł ukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  In _Advances in Neural Information Processing Systems 30_ , pp. 5998–6008. 2017.
  * Voorhees et al. (1999) Ellen M Voorhees et al.  The TREC-8 question answering track report.  In _TREC_ , 1999.
  * Yang & Seo (2020) Sohee Yang and Minjoon Seo.  Is retriever merely an approximator of reader?  _arXiv preprint arXiv:2010.10999_ , 2020.
  * Yang et al. (2019) Wei Yang, Yuqing Xie, Aileen Lin, Xingyu Li, Luchen Tan, Kun Xiong, Ming Li, and Jimmy Lin.  End-to-end open-domain question answering with BERTserini.  In _Proc. NAACL (Demonstrations)_ , 2019.

##  Appendix A Experimental Details

###  A.1 Setting

For NaturalQuestions and TriviaQA we follow the standard open-domain question
answering setting used in Lee et al. (2019); Karpukhin et al. (2020). In this
setting the original development set is used as test set, and 10% of the
training set is used for development purpose. Moreover, for NaturalQuestions,
all questions with answers longer than five tokens are discarded.

For TriviaQA we use the unique human-generated answer to train the reader. In
this dataset part of the answers are in uppercase. We normalize uppercase
answers by converting the first letter in each word to uppercase and remaining
characters to lowercase using the title Python string method.

For NarrativeQA, questions and answers in uppercase are converted to
lowercase.

###  A.2 Training

For every datasets, both the reader and the retriever are fine-tuned with a
dropout rate of 10%. All models at the exception of the large reader are
trained using the ADAM algorithm (Kingma & Ba, 2014) with a constant learning
rate of $10^{-4}$ for the base reader and $5\cdot 10^{-5}$ for the retriever.
The base reader is trained for 10k gradient steps with a batch size of 64. We
train the large reader with the ADAMW algorithm (Loshchilov & Hutter, 2019)
with a peak learning rate of $5\cdot 10^{-5}$ and a linear warmup for 600
gradient steps followed by a linear decrease of the learning rate for 14.4k
gradient steps.

We perform model selection on the validation performance. The retriever is
trained until its performance saturates with a batch size of 64. To monitor
the performance of the retriever during training, we measure the similarity
between the ranking obtained with the reader score, and the ranking of the
retriever. We use different metrics for this: the number of inversions between
the two rankings, the proportion of passages in the retriever top-$k$ that are
also in the reader top-$k$ and the number of passages to obtain all top-$k$
passage of the reader.

During training and at test time, each text input of the encoder is restricted
to be at most 250 token long. For NaturalQuestions and TriviaQA, we use
wikipedia as a knowledge source, thus for each passage there is an associated
article title. Each input is composed of the concatenation of a question,
title and support passage with special tokens question:, title: and context:
added before the question, the title and the text of each passage. In the case
of NarrativeQA, the question and each passage are concatenated to form the
different inputs.

Hyperparameter | Reader-base | Reader-large | Retriever
---|---|---|---
Number of parameters | 220M | 770M | 110M
Number of heads | 12 | 16 | 12
Number of layers | 24 | 48 | 12
Hidden size | 768 | 1024 | 768
Batch size | 64 | 64 | 64
Dropout | 0.1 | 0.1 | 0.1
Learning rate schedule | constant | linear | constant
Peak learning rate | 0.0001 | 0.00005 | 0.00005
Gradient clipping | 1. | 1. | 1.
Table 6: Hyperparameters for retriever and reader training.

###  A.3 Inference

At test time, for TriviaQA and NaturalQuestions we use greedy decoding, and
Beam Search with 3 beams for NarrativeQA.

|  | NaturalQuestions |  | TriviaQA
---|---|---|---|---
Iter. |  | P@20 | P@100 | Dev EM |  | P@20 | P@100 | Dev EM
0 |  | 77.1 | 84.3 | 46.4 |  | 78.2 | 84.7 | 65.0
1 |  | 80.3 | 86.7 | 47.8 |  | 81.4 | 86.4 | 67.1
2 |  | 82.4 | 87.9 | 48.2 |  | 83.5 | 87.4 | 68.1
Table 7: Iterative training starting with documents retrieved with DPR.
Iteration 0 corresponds to the performance of the reader trained on the set of
initial support documents. We report all metrics on the validation set.
Contrary to results reported in Table 1, the reader model was not re-
initialized between each iteration.


# Computing the Barnes $G$-function and the gamma function in the entire
complex plane

Alexey Kuznetsov 111Dept. of Mathematics and Statistics, York University, 4700
Keele Street, Toronto, ON, M3J 1P3, Canada. Email<EMAIL_ADDRESS>Webpage:
kuznetsov.mathstats.yorku.ca

###### Abstract

We present an algorithm for generating approximations for the logarithm of
Barnes $G$-function in the half-plane $\textnormal{Re}(z)\geq 3/2$. These
approximations involve only elementary functions and are easy to implement.
The algorithm is based on a two-point Padé approximation and we use it to
provide two approximations to $\ln(G(z))$, accurate to $3\times 10^{-16}$ and
$3\times 10^{-31}$ in the half-plane $\textnormal{Re}(z)\geq 3/2$; a
reflection formula is then used to compute Barnes $G$-function in the entire
complex plane. A by-product of our algorithm is that it also produces accurate
approximations to the gamma function.

Keywords: Barnes $G$-function, gamma function, Laplace transform, Pade
approximation
2010 Mathematics Subject Classification : Primary 33F05, Secondary 30E10,
33E30, 33B15

##  1 Introduction

Barnes $G$-function is an entire function defined via the Weierstrass product

$G(1+z)=(2\pi)^{z/2}e^{-\frac{1}{2}(z+(1+\gamma)z^{2})}\prod\limits_{k\geq
1}\Big{(}1+\frac{z}{k}\Big{)}^{k}e^{-z+\frac{z^{2}}{2k}},\;\;\;z\in{\mathbb{C}}.$
(1)

Here $\gamma$ is the Euler-Mascheroni constant. This function satisfies

$G(z+1)=\Gamma(z)G(z).$ (2)

and it is known [22] that $G$ is a unique function satisfying (2) subject to
conditions $G(1)=1$ and

$\frac{d^{3}}{dz^{3}}\ln G(z)\geq 0\;\;{\textrm{for }}z>1.$

The systematic study of $G(z)$ (and certain related functions, including a
more general double gamma function $G(z;\tau)$ and a hierarchy of multiple
gamma functions $\Gamma_{n}(z)$) was initiated by Barnes in a series of papers
[2, 3, 4]. The $G$-function was in fact introduced in an earlier paper by
Alexeiewsky [1] and in some older literature it was sometimes called
“Alexeiewsky function”, though by now the name “Barnes $G$-function” seems to
be universally accepted.

Starting in 1980s the Barnes $G$-function and related double gamma functions
started to appear frequently in various areas of Mathematics and Physics.
Shintani [20] used double gamma function to give a new proof of Kronecker’s
limit formula. The Barnes $G$-function was found to have connections to
determinants of Laplacians [19, 21] and it also appeared in the Keating-Snaith
conjecture on moments of the Riemann zeta function [10]. More recently this
function has found many aplications in Probability Theory [11, 16, 23].

Despite this increased interest in the Barnes $G$-function and its apparent
usefullness in many applications, no work was done so far on developing
efficient numerical algorithms for computing $G(z)$ (though several recent
papers studied asymptotic expansions of $G(z)$, see [8, 9, 15, 24]). It is the
purpose of this paper to correct this omission. Our main contribution is an
algorithm that can produce accurate and simple approximations to $\ln(G(z))$
in the half-plane $\textnormal{Re}(z)\geq 3/2$ – we demonstrate this by
providing two such approximations, one valid to $3\times 10^{-16}$ and another
to $3\times 10^{-31}$. With the help of a reflection formula and the
functional equation (2) this allows one to compute $\ln(G(z))$ in the entire
complex plane. Our algorithm also generats accurate approximations to the
gamma function, though this is probably less important as there already exist
excellent approximations to $\Gamma(z)$, such as the Lanczos approximation
[12, 17].

The key part of our algorithm consists in approximating the function

$f(x):=\frac{e^{-x}}{x^{3}}\left[\frac{1}{2}\coth\left(\frac{x}{2}\right)-\frac{1}{x}-\frac{x}{12}\right],\;\;\;x>0,$
(3)

by a finite sum of exponential functions
$\sum_{j=1}^{p}c_{j}e^{-\lambda_{j}x}$, and this is done by finding a two-
point Padé approximation to the Laplace transform of $f(x)$. This two-point
Padé approximation method is of independent interest, as it can be used to
provide approximations (by sums of exponential functions) to any function $g$
for which derivatives $g^{(n)}(0)$ and “moments”
$\int_{0}^{\infty}x^{n}g(x){\textnormal{d}}x$ are known explicitly or can be
computed to sufficiently high precision.

The paper is organized as follows. In Section 2 we explain how approximations
to $f(x)$ by a finite sum of exponential functions lead to approximations to
$\ln(\Gamma(z))$ and $\ln(G(z))$ in the half-plane $\textnormal{Re}(z)\geq
3/2$ and how one can extend these approximations to the entire complex plane
by using functional equations and reflection formulas. Section 3 describes the
two-point Padé approximation algorithm for finding exponential sum
approximations and in Section 4 we present numerical results that illustrate
the performance of this algorithm. In Appendix A we describe an algorithm for
computing the dilogarithm function (this is needed for implementing the
reflection formula for the Barnes $G$-function), and in Appendix B we present
the numerical values of the coefficients $c_{j}$ and $\lambda_{j}$.

##  2 Approximatng logarithms of $\Gamma(z)$ and $G(z)$

The starting point for our investigation is the following proposition, which
reduces the computation of the gamma function and the Barnes $G$-function to
the computation of

$F(z):=\int\limits_{0}^{\infty}e^{-zx}xf(x){\textnormal{d}}x$ (4)

and its derivative $F^{\prime}(z)$. In what follows, $A$ denotes the the
Glaisher-Kinkelin constant and $\ln(z)$ denotes the principal branch of the
logarithm function.

######  Proposition 1.

For $\textnormal{Re}(z)>0$

$\displaystyle\ln(\Gamma(z))$
$\displaystyle=\left(z-\frac{1}{2}\right)\ln(z)-z+\frac{1}{2}\ln(2\pi)+\frac{1}{12z}-F^{\prime}(z-1),$
(5) $\displaystyle\ln(G(z))$
$\displaystyle=\left(\frac{z^{2}}{2}-z+\frac{5}{12}\right)\ln(z)-\frac{3}{4}z^{2}+\frac{1}{2}\ln(2\pi)(z-1)+z$
(6)
$\displaystyle+\frac{1}{12}-\ln(A)-\frac{1}{12z}+F(z-1)-(z-1)F^{\prime}(z-1).$

###### Proof.

Formula (5) is obtained from the well-known Binet’s integral representation

$\ln(\Gamma(z))=\left(z-\frac{1}{2}\right)\ln(z)-z+\frac{1}{2}\ln(2\pi)+\int_{0}^{\infty}\frac{e^{-zx}}{x}\left[\frac{1}{e^{x}-1}-\frac{1}{x}+\frac{1}{2}\right]{\textnormal{d}}x,$

by noting that $f$ can be expressed in an equivalent form

$f(x)=\frac{e^{-x}}{x^{3}}\left[\frac{1}{e^{x}-1}-\frac{1}{x}+\frac{1}{2}-\frac{x}{12}\right].$
(7)

Next, by setting $N=1$ in [9][Equation (12)] we obtain

$\ln(G(z))=(z-1)\ln(\Gamma(z))-\frac{1}{2}\left(z^{2}-z+\frac{1}{6}\right)\ln(z)+\frac{z^{2}}{4}-\ln(A)+F(z-1),$

and this result, together with (5), implies (6).
$\sqcap\kern-8.0pt\hbox{$\sqcup$}$

Next we describe our approximation algorithm.

### Approximations in the domain $\textnormal{Re}(z)\geq 3/2$

Let $p$ be a positive integer. Assume we have found coefficients
$\\{c_{j}\\}_{1\leq j\leq p}$ and $\\{\lambda_{j}\\}_{1\leq j\leq p}$ such
that $\textnormal{Re}(\lambda_{j})>0$ and

$\phi(x):=\sum\limits_{j=1}^{p}c_{j}e^{-\lambda_{j}x}$ (8)

is a “good” approximation to $f(x)$ on $x\in(0,\infty)$, in the sense that the
numbers

$\displaystyle\epsilon_{1}$
$\displaystyle:=2\sup\limits_{x>0}\Big{|}x^{2}(f(x)-\phi(x))\Big{|},$ (9)
$\displaystyle\epsilon_{2}$
$\displaystyle:=2\sup\limits_{x>0}\Big{|}3x(f(x)-\phi(x))+x^{2}(f^{\prime}(x)-\phi^{\prime}(x))\Big{|},$
(10)

are small. For $\textnormal{Re}(z)>0$ we define a rational function

$\Phi(z):=\sum\limits_{j=1}^{p}\frac{c_{j}}{(z+\lambda_{j})^{2}}=\int_{0}^{\infty}e^{-zx}x\phi(x){\textnormal{d}}x$
(11)

and we define approximations to $\ln(\Gamma(z))$ and $\ln(G(z))$ in the form

$\displaystyle\ln(\widehat{\Gamma}(z))$
$\displaystyle=\left(z-\frac{1}{2}\right)\ln(z)-z+\frac{1}{2}\ln(2\pi)+\frac{1}{12z}-\Phi^{\prime}(z-1),$
(12) $\displaystyle\ln(\widehat{G}(z))$
$\displaystyle=\left(\frac{z^{2}}{2}-z+\frac{5}{12}\right)\ln(z)-\frac{3}{4}z^{2}+\frac{1}{2}\ln(2\pi)(z-1)+z$
(13)
$\displaystyle+\frac{1}{12}-\ln(A)-\frac{1}{12z}+\Phi(z-1)-(z-1)\Phi^{\prime}(z-1).$

######  Proposition 2.

In the half-plane $\textnormal{Re}(z)\geq 3/2$ we have

$|\ln(\Gamma(z))-\ln(\widehat{\Gamma}(z))|\leq\epsilon_{1}$ (14)

and

$|\ln(G(z))-\ln(\widehat{G}(z))|\leq\epsilon_{2}.$ (15)

###### Proof.

From (5) and (12) we see that

$\ln(\Gamma(z))-\ln(\widehat{\Gamma}(z))=\Phi^{\prime}(z-1)-F^{\prime}(z-1).$

We set $w=z-1$ and estimate for $\textnormal{Re}(w)\geq 1/2$

$|F^{\prime}(w)-\Phi^{\prime}(w)|=\Big{|}\int_{0}^{\infty}e^{-wx}x^{2}(f(x)-\phi(x)){\textnormal{d}}x\Big{|}\leq\int_{0}^{\infty}e^{-x/2}x^{2}|f(x)-\phi(x)|{\textnormal{d}}x\leq\frac{\epsilon_{1}}{2}\int_{0}^{\infty}e^{-x/2}{\textnormal{d}}x=\epsilon_{1},$

which proves (14). Similarly, from (6) and (13) we see that

$\ln(G(z))-\ln(\widehat{G}(z))=F(z-1)-(z-1)F^{\prime}(z-1)-\Phi(z-1)+(z-1)\Phi^{\prime}(z-1).$

We perform integration by parts

$wF^{\prime}(w)=-w\int_{0}^{\infty}e^{-wx}x^{2}f(x){\textnormal{d}}x=-\int_{0}^{\infty}e^{-wx}(2xf(x)+x^{2}f^{\prime}(x)){\textnormal{d}}x$

and estimate for $\textnormal{Re}(w)\geq 1/2$

$\displaystyle|F(w)-wF^{\prime}(w)-\Phi(w)+w\Phi^{\prime}(w)|$
$\displaystyle=\Big{|}\int_{0}^{\infty}e^{-wx}(3x(f(x)-\phi(x))+x^{2}(f^{\prime}(x)-\phi^{\prime}(x))){\textnormal{d}}x\Big{|}$
$\displaystyle\leq\int_{0}^{\infty}e^{-x/2}|3x(f(x)-\phi(x))+x^{2}(f^{\prime}(x)-\phi^{\prime}(x))|{\textnormal{d}}x$
$\displaystyle\leq\frac{\epsilon_{2}}{2}\int_{0}^{\infty}e^{-x/2}{\textnormal{d}}x=\epsilon_{2}.$

$\sqcap\kern-8.0pt\hbox{$\sqcup$}$

### Extending the domain to $\textnormal{Re}(z)\geq 1/2$

In the vertical strip $1/2\leq\textnormal{Re}(z)<3/2$ we compute
$\ln(\Gamma(z))$ via functional equation

$\ln(\Gamma(z))=\ln(\Gamma(z+1))-\ln(z),$

and then we compute $\ln(G(z))$ via

$\ln(G(z))=\ln(G(z+1))-\ln(\Gamma(z)).$

### Approximations in the domain $\textnormal{Re}(z)<1/2$: applying reflection
formulas

It is well-known that Gamma function satisfies the reflection identity

$\Gamma(z)\Gamma(1-z)=\frac{\pi}{\sin(\pi z)}.$ (16)

In the previous sections we described an algorithm for computing
$\ln(\Gamma(z))$ for $\textnormal{Re}(z)\geq 1/2$ and the above identity gives
a simple algorithm for computing $\Gamma(z)$ for $\textnormal{Re}(z)<1/2$.
However, some additional work is needed if we want to compute $\ln(\Gamma(z))$
for $\textnormal{Re}(z)<1/2$ while preserving analyticity of this function in
${\mathbb{C}}\setminus(-\infty,0]$. First of all, we rewrite the above
reflection identity in the form

$\Gamma(z)=\frac{1}{\Gamma(1-z)}\times\frac{-2\pi{\textnormal{i}}e^{\pi{\textnormal{i}}z}}{1-e^{2\pi{\textnormal{i}}z}}.$

We assume that $\textnormal{Im}(z)>0$ and take logarithms of both sides of the
above equation and obtain

$\ln(\Gamma(z))=-\ln(\Gamma(1-z))+\ln(2\pi)-\frac{\pi{\textnormal{i}}}{2}+\pi{\textnormal{i}}z-\ln(1-e^{2\pi{\textnormal{i}}z}).$
(17)

One can check that we chose the correct branch of logarithm by noting that the
above equation is valid if $\textnormal{Re}(z)=1/2$ and $\textnormal{Im}(z)\to
0^{+}$. The right-hand side of (17) is analytic in $\textnormal{Im}(z)>0$,
thus we can use this identity to compute $\ln(\Gamma(z))$ for
$\textnormal{Re}(z)<1/2$ and $\textnormal{Im}(z)>0$. This identity is also
used for $0<z<1/2$ (again, we need to use the principal branch of the
logarithm). When $\textnormal{Re}(z)<1/2$ and $\textnormal{Im}(z)<0$ we use
the above method coupled with the identity
$\overline{\ln(\Gamma(z))}=\ln(\Gamma(\bar{z}))$.

Computing $\ln(G(z))$ in the half-plane $\textnormal{Re}(z)<1/2$ follows a
similar procedure. In this case the relevant reflection formula is

$\frac{G(1+z)}{G(1-z)}=(2\pi)^{z}\exp\Big{(}-\int_{0}^{z}\pi t\cot(\pi
t){\textnormal{d}}t\Big{)}.$

Using formulas (4.6), (8.43) in [13] the integral in the right-hand side of
the above identity can be expressed in terms of the dilogarithm function

${\textrm{Li}}_{2}(z):=-\int_{0}^{z}\ln(1-t)\frac{{\textnormal{d}}t}{t}=\sum\limits_{k\geq
1}\frac{z^{2}}{k^{2}},\;\;\;|z|\leq 1,$ (18)

so that changing the variable $z\mapsto z-1$ we would obtain a reflection
formula in the form

$\displaystyle\ln(G(z))$
$\displaystyle=\ln(G(2-z))+(z-1)\ln(2\pi)+\frac{\pi{\textnormal{i}}}{2}(z^{2}-2z+5/6)$
(19)
$\displaystyle-(z-1)\ln(1-e^{2\pi{\textnormal{i}}z})-\frac{1}{2\pi{\textnormal{i}}}{\textrm{Li}}_{2}(e^{2\pi{\textnormal{i}}z}),\;\;\;\textnormal{Im}(z)>0.$

Again, one can check that we have the correct branch of the logarithm by
taking $\textnormal{Re}(z)=1$ and $\textnormal{Im}(z)\to 0^{+}$ and using the
fact that $G(1)=1$ and ${\textrm{Li}}_{2}(1)=\pi^{2}/6$. We note that the
right-hand side of (19) is analytic in the half-plane $\textnormal{Im}(z)>0$,
thus we can use (19) to compute $\ln(G(z))$ for $\textnormal{Re}(z)<1/2$ and
$\textnormal{Im}(z)>0$. An efficient algorithm for computing the dilogarithm
function ${\textrm{Li}}_{2}(z)$ for $|z|\leq 1$ is presented in the Appendix
A. We also use formula (19) when $0<z<1/2$. When $\textnormal{Re}(z)<1/2$ and
$\textnormal{Im}(z)<0$ we use the above method coupled with the identity
$\overline{\ln(G(z))}=\ln(G(\bar{z}))$.

This concludes our description of the method for computing $\ln(\Gamma(z))$
and $\ln(G(z))$ in ${\mathbb{C}}\setminus(-\infty,0]$. The only outstanding
question is how to approximate the function $f$ defined in (3) by a finite sum
of exponential functions of the form (8) so that the errors $\epsilon_{1}$ and
$\epsilon_{2}$ in (9) and (10) are small. We propose a solution to this
problem in the next section.

##  3 Approximating a function by a finite sum of exponential functions

There exist a number of methods for approximating a given function by a finite
sum of exponential functions. In [18] the problem of finding the best
approximation on a finite interval in the Chebyshev sense is considered. A
method based on Hankel matrices is introduced in [5] and a Poisson summation
approach is presented in [6]. See also [7] for results on approximating $1/x$
on $[1,\infty)$. In this paper we propose a different approach for
approximating functions by finite sums of exponential functions: this method
is based on a two-point Padé approximation and it is particularly well-suited
for the function $f$ given in (3).

First we show that for the function $f$ defined by (3) one can compute
explicitly (in terms of values of the Riemann zeta function $\zeta(s)$ at
positive integers and some other well-known constants) both the derivatives at
zero and the moments.

######  Lemma 1.

  * (i)

For $n\geq 0$

$f^{(n)}(0)=(-1)^{n}n!\sum\limits_{k=0}^{[n/2]}\frac{B_{2k+4}}{(n-2k)!(2k+4)!},$
(20)

where $B_{n}$ are Bernoulli numbers.

  * (ii)

Let $\mu_{n}:=\int_{0}^{\infty}x^{n}f(x){\textnormal{d}}x$. Then

$\displaystyle\mu_{0}$
$\displaystyle=\frac{1}{2}\zeta^{\prime}(-2)+\frac{1}{72}=-\frac{\zeta(3)}{8\pi^{2}}+\frac{1}{72},$
$\displaystyle\mu_{1}$
$\displaystyle=-\zeta^{\prime}(-1)-\frac{1}{6}=\ln(A)-\frac{1}{4},$
$\displaystyle\mu_{2}$
$\displaystyle=\zeta^{\prime}(0)+\frac{11}{12}=-\ln(2\pi)/2+\frac{11}{12},$
$\displaystyle\mu_{3}$ $\displaystyle=\lim\limits_{s\to
1}\Big{(}\zeta(s)-\frac{1}{s-1}\Big{)}-\frac{7}{12}=\gamma-7/12,$

and for $n\geq 4$

$\mu_{n}=(n-3)!\Big{[}\zeta(n-2)-\frac{1}{n-3}-\frac{1}{2}-\frac{n-2}{12}\Big{]}.$

Here $\gamma$ is the Euler-Mascheroni constant and $A$ is the Glaisher-
Kinkelin constant.

###### Proof.

Formula (20) follows from the Taylor series expansion of the hyperbolic
cotangent function:

$\frac{1}{2}\coth\left(\frac{x}{2}\right)=\sum\limits_{n=0}^{\infty}\frac{B_{2n}}{(2n)!}x^{2n-1}.$

Formulas for moments $\mu_{n}$ follow from the identity

$\int_{0}^{\infty}x^{s+2}f(x){\textnormal{d}}x=\Gamma(s)\Big{[}\zeta(s)-\frac{1}{s-1}-\frac{1}{2}-\frac{s}{12}\Big{]},\;\;\;\textnormal{Re}(s)>-3,$
(21)

which can be obtained from (7) and integral representations

$\displaystyle\int_{0}^{\infty}\frac{e^{-x}x^{s-1}}{e^{x}-1}{\textnormal{d}}x=\Gamma(s)(\zeta(s)-1),\;\;\;\textnormal{Re}(s)>1,$
$\displaystyle\int_{0}^{\infty}e^{-x}x^{s-1}{\textnormal{d}}x=\Gamma(s),\;\;\;\textnormal{Re}(s)>0.$

The above two identities can be used to deduce (21) in the domain
$\textnormal{Re}(s)>1$ and extension of (21) to the half-plane
$\textnormal{Re}(s)>-3$ is done by analytic continuation (note that both sides
of (21) are well-defined and analytic in $\textnormal{Re}(s)>-3$). Now we
compute

$\displaystyle\mu_{0}=\int_{0}^{\infty}f(x){\textnormal{d}}x$
$\displaystyle=\lim\limits_{s\to-2}\Gamma(s)\Big{[}\zeta(s)-\frac{1}{s-1}-\frac{1}{2}-\frac{s}{12}\Big{]}$
$\displaystyle=\lim\limits_{s\to-2}(s-2)\Gamma(s)\times\frac{\zeta(s)-\frac{1}{s-1}-\frac{1}{2}-\frac{s}{12}}{s-2}$
$\displaystyle=\frac{1}{2}\frac{{\textnormal{d}}}{{\textnormal{d}}s}\Big{[}\zeta(s)-\frac{1}{s-1}-\frac{1}{2}-\frac{s}{12}\Big{]}\bigg{|}_{s=-2}=\frac{1}{2}\zeta^{\prime}(-2)+\frac{1}{72},$

and the second expression for $\mu_{0}$ is obtained from the reflection
formula: $\zeta^{\prime}(-2)=-\zeta(3)/(4\pi^{2})$. The proof of formulas for
$\mu_{1}$, $\mu_{2}$ and $\mu_{3}$ is identical and we leave all details to
the reader. Formulas for $\mu_{n}$ with $n\geq 4$ follow directly from (21).
$\sqcap\kern-8.0pt\hbox{$\sqcup$}$

Now we describe our method for finding approximations to $f(x)$ in the form
$\phi(x)=\sum\limits_{j=1}^{p}c_{j}e^{-\lambda_{j}x}$. We fix two non-negative
integers $N_{1}$ and $N_{2}$ such that $N_{1}+N_{2}=2p$ and we aim to find
coefficients $\\{c_{j}\\}_{1\leq j\leq p}$ and $\\{\lambda_{j}\\}_{1\leq j\leq
p}$ such that the functions $f$ and $\phi$ have the same first $N_{1}$ moments
and the same first $N_{2}$ terms in Taylor series expansion at zero. In other
words, we want to enforce the following conditions:

$\displaystyle\int_{0}^{\infty}x^{n}\phi(x){\textnormal{d}}x=n!\sum\limits_{j=1}^{p}c_{j}\lambda_{j}^{-n-1}=\mu_{n},\;\;\;n=0,1,2,\dots,N_{1}-1,$
(22)
$\displaystyle\phi^{(n)}(0)=(-1)^{n}\sum\limits_{j=1}^{p}c_{j}\lambda_{j}^{n}=f^{(n)}(0),\;\;\;n=0,1,2,\dots,N_{2}-1.$
(23)

As we will show next, this problem is equivalent to a two-point Padé
approximation, with the two points being $z=0$ and $z=\infty$. Indeed, let us
consider the Laplace transform

$\Psi(z):=\int_{0}^{\infty}e^{-zx}\phi(x){\textnormal{d}}x=\sum\limits_{j=1}^{p}\frac{c_{j}}{z+\lambda_{j}}.$
(24)

From the above partial fraction decomposition it is easy to obtain the
following series representations

$\displaystyle\Psi(z)$ $\displaystyle=\sum\limits_{n\geq
0}z^{-n-1}\Big{[}(-1)^{n}\sum_{j=1}^{p}c_{j}\lambda_{j}^{n}\Big{]},$
$\displaystyle\Psi(z)$ $\displaystyle=\sum\limits_{n\geq
0}z^{n}\Big{[}(-1)^{n}\sum_{j=1}^{p}c_{j}\lambda_{j}^{-n-1}\Big{]},$

where the first one converges for all $z$ large enough and the second for all
$z$ in a neighbourhood of zero. In light of the above equations we see that
the problem of finding a function
$\phi(x)=\sum\limits_{j=1}^{p}c_{j}e^{-\lambda_{j}x}$ that satisfies
conditions (22) and (23) is equivalent to the following problem: find the
coefficients $\\{a_{j}\\}_{1\leq j\leq p}$ and $\\{b_{j}\\}_{1\leq j\leq p}$
of the rational function

$\Psi(z)=\frac{P(z)}{Q(z)}=\frac{a_{0}+a_{1}z+\dots+a_{p-1}z^{p-1}}{b_{0}+b_{1}z+\dots+b_{p-1}z^{p-1}+z^{p}}$
(25)

such that

$\displaystyle\Psi(z)$
$\displaystyle=\beta_{0}+\beta_{1}z+\beta_{2}z^{2}+\dots+\beta_{N_{1}-1}z^{N_{1}-1}+O(z^{N_{1}}),\;\;\;z\to
0,$ (26) $\displaystyle\Psi(z)$
$\displaystyle=\alpha_{0}z^{-1}+\alpha_{1}z^{-2}+\dots\alpha_{N_{2}-1}z^{-N_{2}}+O(z^{-N_{2}-1}),\;\;\;z\to\infty,$
(27)

where we denoted

$\alpha_{n}:=f^{(n)}(0),\;\;\;\beta_{n}:=\frac{(-1)^{n}}{n!}\mu_{n}.$ (28)

This is precisely the two-point Padé approximation problem. As in the case of
the classical (one-point) Padé approximation, the above problem can be reduced
to solving a system of linear equations. First we rewrite equation (26) in the
form

$\displaystyle a_{0}+a_{1}z+\dots+a_{p-1}z^{p-1}$
$\displaystyle=\Big{(}b_{0}+b_{1}z+\dots+b_{p-1}z^{p-1}+z^{p}\Big{)}$
$\displaystyle\times\Big{(}\beta_{0}+\beta_{1}z+\beta_{2}z^{2}+\dots+\beta_{N_{1}-1}z^{N_{1}-1}+O(z^{N_{1}})\Big{)}.$

Comparing the coefficients in front of $z^{k}$ with $k=0,1,\dots,N_{1}-1$ we
obtain equations

$a_{k}{\mathbf{1}}_{\\{k\leq p-1\\}}=\beta_{k-p}{\mathbf{1}}_{\\{p\leq
k\\}}+\sum_{j=0}^{p-1}b_{j}\beta_{k-j}{\mathbf{1}}_{\\{0\leq k-j\leq
N_{1}-1\\}}.$ (29)

Next we introduce a variable $u=1/z$ and rewrite (27) in the form

$\displaystyle a_{p-1}+a_{p-2}u+\dots+a_{0}u^{p-1}$
$\displaystyle=\Big{(}1+b_{p-1}u+\dots+b_{0}u^{p}\Big{)}$
$\displaystyle\times\Big{(}\alpha_{0}+\alpha_{1}u+\dots+\alpha_{N_{2}-1}u^{N_{2}-1}+O(u^{N_{2}})\Big{)}.$

Comparing the coefficients in front of $u^{k}$ with $k=0,1,\dots,N_{2}-1$ we
obtain equations

$a_{p-1-k}{\mathbf{1}}_{\\{k\leq
p-1\\}}=\alpha_{k}+\sum\limits_{j=0}^{p-1}b_{j}\alpha_{k+j-p}{\mathbf{1}}_{\\{0\leq
k+j-p\leq N_{2}-1\\}}.$ (30)

Now we have $N_{1}$ linear equations of the form (29) and $N_{2}$ linear
equations of the form (30) – together these give us $N_{1}+N_{2}=2p$ equations
for $2p$ unknowns $\\{a_{j},b_{j}\\}_{0\leq j\leq p-1}$. To rewrite these
equations in a matrix form (that is more suitable for numerical computation),
we introduce a $2p\times 2p$ matrix ${\bf A}=\\{A_{i,j}\\}_{1\leq i,j\leq 2p}$
and a column vector $\vec{\bf y}\in{\mathbb{R}}^{2p}$ by specifying their
elements as follows: for $k=0,1,\dots,N_{1}-1$ and $j=0,1,\dots,p-1$

$\displaystyle A_{k+1,k+1}$ $\displaystyle={\mathbf{1}}_{\\{k\leq p-1\\}},$
$\displaystyle A_{k+1,p+1+j}$
$\displaystyle=-\beta_{k-j}{\mathbf{1}}_{\\{0\leq k-j\leq N_{1}-1\\}},$
$\displaystyle y_{k+1}$ $\displaystyle=\beta_{k-p}{\mathbf{1}}_{\\{p\leq
k\\}},$

and for $l=0,1,\dots,N_{2}-1$ and $j=0,1,\dots,p-1$

$\displaystyle A_{N_{1}+l+1,p-l}$ $\displaystyle={\mathbf{1}}_{\\{l\leq
p-1\\}},$ $\displaystyle A_{N_{1}+l+1,p+1+j}$
$\displaystyle=-\alpha_{l+j-p}{\mathbf{1}}_{\\{0\leq l+j-p\leq N_{2}-1\\}},$
$\displaystyle y_{N_{1}+l+1}$ $\displaystyle=\alpha_{l}.$

Now the two sets of linear equations (29) and (30) can be represented as a
single matrix equation ${\bf A}\vec{\bf x}=\vec{\bf y}$, where a column vector
$\vec{\bf x}\in{\mathbb{R}}^{2p}$ is defined as

$\vec{\bf x}:=[a_{0},a_{1},\dots,a_{p-1},b_{0},b_{1},\dots,b_{p-1}]^{T}.$ (31)

After we solve equation ${\bf A}\vec{\bf x}=\vec{\bf y}$ and find $\vec{\bf
x}$, we compute $a_{j}=x_{j+1}$ and $b_{j}=x_{p+j+1}$ for $j=0,1,\dots,p-1$,
and this gives us the desired rational function $\Psi$ in the form (25).
Performing the partial fraction decomposition (and assuming that all the roots
of the denominator $Q(z)$ are simple) we would obtain the expression of $\phi$
in the form (24) and this would give us the desired coefficients
$\\{c_{j}\\}_{1\leq j\leq p}$ and $\\{\lambda_{j}\\}_{1\leq j\leq p}$ of
$\phi(x)=\sum\limits_{j=1}^{p}c_{j}e^{-\lambda_{j}x}$. We can check the
accuracy of computation of these coefficients via (22) and (23).

Remark 1: Note that the functions $\Phi$ and $\Psi$ defined in (11) and (24)
satisfy $\Phi(z)=-\Psi^{\prime}(z)$. Thus the rational function $\Phi$ has the
same $N_{1}-1$ first terms of Taylor expansion at zero and the same first
$N_{2}+1$ asymptotic terms as $z\to+\infty$ as the function $F(z)$ . Therefore
our approximations $\ln(\widehat{\Gamma}(z))$ and $\ln(\widehat{G}(z))$, which
were introduced in (14) and (15), satisfy the following properties: as $z\to
1^{+}$

$\displaystyle\ln(\Gamma(z))-\ln(\widehat{\Gamma}(z))=O((z-1)^{N_{1}-1}),$
$\displaystyle\ln(G(z))-\ln(\widehat{G}(z))=O((z-1)^{N_{1}-1}),$

and as $z\to\infty$

$\displaystyle\ln(\Gamma(z))-\ln(\widehat{\Gamma}(z))=O(z^{-N_{2}-2}),$
$\displaystyle\ln(G(z))-\ln(\widehat{G}(z))=O(z^{-N_{2}-2}).$

The first of these is probably not very informative, as we only use
approximations $\ln(\widehat{\Gamma}(z))$ and $\ln(\widehat{G}(z))$ in the
half-plane $\textnormal{Re}(z)\geq 3/2$ (where $z$ is bounded away from $1$),
but the second result is useful as it tells us that our approximations in the
half-plane $\textnormal{Re}(z)\geq 3/2$ become more accurate as $|z|$
increases.

To concluse this section, we would like to point out that there is no
theoretical justification for the above algorithm. One can not be certain that
the system ${\bf A}\vec{\bf x}=\vec{\bf y}$ has a solution or that all roots
of polynomial $Q(z)$ in (25) are simple (though in all our numerical
experiments this was always the case). Even if the algorithm does produce the
numbers $\\{c_{j}\\}_{1\leq j\leq p}$ and $\\{\lambda_{j}\\}_{1\leq j\leq p}$,
there is no guarantee that this would result in a good approximation
$\phi(x)=\sum\limits_{j=1}^{p}c_{j}e^{-\lambda_{j}x}$ to $f(x)$ (in the sense
of errors $\epsilon_{1}$ and $\epsilon_{2}$ in (9) and (10) being small).
Existing results on convergence of Padé approximations to Stieltjes functions
are not applicable here since the function $f(x)$ is not completely monotone
(thus its Laplace transform is not a Stieltjes function). Lacking any
theoretical justification for the above algorithm, in the next section we will
try to convince the reader that it works well by presenting numerical
evidence.

##  4 Numerical results

First we would like to make several remarks about implementation of the
algorithm presented in Section 3. The algorithm includes two steps that are
ill-conditioned and that require high precision arithmetic: (i) solving the
system of linear equations ${\bf A}\vec{\bf x}=\vec{\bf y}$ and (ii) finding
the partial fraction decomposition of the rational function $\Psi(z)$ (which
is equivalent to finding roots of the polynomial $Q(z)$ in (25)). Because of
this issue, we perform all computations with the precision of 500 digits (we
use David Bailey’s Fortran module MPFUN). The inputs of the algorithm are the
constants $f^{(n)}(0)$ and $\mu_{n}$ given in Lemma 1: we pre-compute them
also to the precision of 500 digits. We used existing tables of Bernoulli
numbers (which also give us the values of $\zeta(2n)$) and we computed the
values of the $\zeta(2n+1)$ via Euler-Maclaurin summation.

(a) The graph of $f(x)$

(b) The graph of $f(x)-\sum\limits_{j=1}^{3}c_{j}e^{-\lambda_{j}x}$

Figure 1: Here $\\{\lambda_{j},c_{j}\\}_{1\leq j\leq 3}$ are obtained by a
two-point Padé approximation algorithm with $p=N_{1}=N_{2}=3$.

First we do a test run of our approximation algorithm with $p=N_{1}=N_{2}=3$.
The two-point Padé approximation algorithm produces the following values of
$\\{\lambda_{i}\\}_{1\leq i\leq 3}$ and $\\{c_{i}\\}_{1\leq i\leq 3}$:

$\displaystyle\lambda_{1}$ $\displaystyle=1.119578657398,$
$\displaystyle\lambda_{2}$
$\displaystyle=1.348997960953+{\textnormal{i}}\times 0.462842303653,$
$\displaystyle c_{1}$ $\displaystyle=-1.501034759773\times 10^{-3},$
$\displaystyle c_{2}$ $\displaystyle=5.607293544244\times
10^{-5}-{\textnormal{i}}\times 1.516208003305\times 10^{-4},$
$\displaystyle\lambda_{3}$
$\displaystyle={\overline{\lambda_{2}}},\;\;\;c_{3}={\overline{c_{2}}}.$

One can check (numerically) that these numbers indeed satisfy conditions (22)
and (23), thus we have found a function
$\phi(x)=\sum\limits_{j=1}^{3}c_{j}e^{-\lambda_{j}x}$ which has the same first
three coefficients of Taylor expansion at zero as $f(x)$ and the same first
three moments as $f(x)$. We present the graphs of $f(x)$ and of the difference
$f(x)-\phi(x)$ on Figure 1. The results already look promising: we can
approximate $f(x)$ to accuracy of $2.2\times 10^{-8}$ with a sum of just three
exponential functions.

| $N_{1}=0$ | $N_{1}=1$ | $N_{1}=2$ | $N_{1}=3$ | $N_{1}=4$ | $N_{1}=5$
---|---|---|---|---|---|---
$N_{2}=0$ |  |  | $4.3\times 10^{-5}$ |  | $4.0\times 10^{-6}$ |
$N_{2}=1$ |  | $1.2\times 10^{-5}$ |  | $6.1\times 10^{-7}$ |  | $9.7\times 10^{-9}$
$N_{2}=2$ | $1.6\times 10^{-5}$ |  | $4.5\times 10^{-7}$ |  | $1.5\times 10^{-8}$ |
$N_{2}=3$ |  | $7.2\times 10^{-7}$ |  | $2.2\times 10^{-8}$ |  | $7.4\times 10^{-9}$
$N_{2}=4$ | $2.3\times 10^{-6}$ |  | $3.7\times 10^{-8}$ |  | $3.6\times 10^{-9}$ |
$N_{2}=5$ |  | $1.1\times 10^{-7}$ |  | $4.2\times 10^{-9}$ |  | $3.5\times 10^{-10}$
Table 1: The error or approximation $\sup\limits_{x\geq
0}|f(x)-\sum\limits_{j=1}^{p}c_{j}e^{-\lambda_{j}x}|$ for various values of
$N_{1}$ and $N_{2}$.

Next we performed a number of similar experiments with all values of $N_{1}$
and $N_{2}$ between zero and five (restricted by the requirement that
$p=(N_{1}+N_{2})/2$ must be an integer). In each case we computed numerically
the $L_{\infty}$ norm

$\sup\limits_{x\geq 0}|f(x)-\sum\limits_{j=1}^{p}c_{j}e^{-\lambda_{j}x}|.$

The results of these computations can be seen in Table 1. We see that the
accuracy of approximation improves with increasing $p$ or $N_{1}$ or $N_{2}$.

(a) $N_{1}=22$, $N_{2}=8$ and $p=15$

(b) $N_{1}=55$, $N_{2}=35$ and $p=45$

Figure 2: The graphs of the functions $\eta_{1}(x)$ (black line) and
$\eta_{2}(x)$ (gray line with circles).

Our next goal was to find values of $N_{1}$ and $N_{2}$ which generate
approximations $\phi(x)=\sum\limits_{j=1}^{p}c_{j}e^{-\lambda_{j}x}$ to $f(x)$
such that the errors $\epsilon_{1}$ and $\epsilon_{2}$ defined in (9) and (10)
are small. We aimed to keep $p$ as small as possible (so that the number of
terms required in the computation of $\Phi(z)$ in (11) is small) and to have
the errors $\epsilon_{1}$ and $\epsilon_{2}$ in the range
$(10^{-16},10^{-15})$, so that this would be adequate for double precision
implementation of approximations to $\ln(\Gamma(z))$ and $\ln(G(z))$. In
another experiment we aimed for the range $(10^{-31},10^{-30})$, that would be
suitable for quadruple precision implementation of $\ln(\Gamma(z))$ and
$\ln(G(z))$. By trial and error we found that taking $p=15$, $N_{1}=22$ and
$N_{2}=8$ produces an approximation $\phi(x)$ with $\epsilon_{1}<10^{-16}$ and
$\epsilon_{2}<3\times 10^{-16}$. These upper bonds are not rigorous, rather
they were established numerically by examining the graphs of

$\displaystyle\eta_{1}(x)$ $\displaystyle:=x^{2}(f(x)-\phi(x)),$
$\displaystyle\eta_{2}(x)$
$\displaystyle:=3x(f(x)-\phi(x))+x^{2}(f^{\prime}(x)-\phi^{\prime}(x)).$

Note that, according to definitions (9) and (10), we have

$\epsilon_{i}=2\sup\limits_{x>0}|\eta_{i}(x)|.$

On Figure 2a the graph of $\eta_{1}(x)$ is shown as a solid black line and the
graph of $\eta_{2}(x)$ as a gray line marked with circles. Figure 2a provides
convincing evidence that $|\eta_{1}(x)|<5\times 10^{-17}$ for all values of
$x$, which implies $\epsilon_{1}<10^{-16}$. The estimate $\epsilon_{2}<3\times
10^{-16}$ is established in the same way by inspecting the graph of
$\eta_{2}(x)$. The coefficients $\\{\lambda_{j}\\}_{1\leq j\leq 15}$ and
$\\{c_{j}\\}_{1\leq j\leq 15}$ obtained by the two-point Padé algorithm in
this case $N_{1}=22$ and $N_{2}=8$ are presented in Appendix B and can also be
downloaded as a list in a text file from the author’s webpage. Using these
coefficients we have implemented in Matlab the approximation algorithm for
$\ln(G(z))$ that was described in Section 2 – this Matlab program can also be
downloaded from the author’s webpage. We used this Matlab program to create
the surface plots of real and imaginary part of $\ln(G(z))$ that are shown on
Figure 3.

Similarly, when $p=45$, $N_{1}=55$ and $N_{2}=35$ we present the graphs of
$\eta_{1}(x)$ and $\eta_{2}(x)$ in Figure 2b and from these graphs we estimate
that $\epsilon_{1}<10^{-31}$ and $\epsilon_{2}<3\times 10^{-31}$. The
coefficients $\\{\lambda_{j}\\}_{1\leq j\leq 45}$ and $\\{c_{j}\\}_{1\leq
j\leq 45}$ corresponding to the case $N_{1}=55$ and $N_{2}=35$ can be found in
the arXiv version of the paper and can also be downloaded from the author’s
webpage. These coefficients may be used to implement quadruple precision
approximations to the logarithm of the Barnes $G$-function, giving more than
30 correct digits in the half-plane $\textnormal{Re}(z)\geq 3/2$.

(a) $\textnormal{Re}\ln(G(z))$

(b) $\textnormal{Im}\ln(G(z))$

Figure 3: The surface plots of real and imaginary parts of $\ln(G(z))$.

## Acknowledgements

Research was supported by the Natural Sciences and Engineering Research
Council of Canada. The author would like to thank two anonymous referees for
careful reading of the paper and for providing constructive comments and
suggestions.

## References

  * [1] W.P. Alexeiewsky.  Uber eine classe von funktionen, die der gammafunktion analog sind.  Leipzig Weidmannsche Buchhandlung, 46:268–275, 1894.
  * [2] E.W. Barnes.  The genesis of the double gamma function.  Proc. London Math. Soc., 31:358–381., 1899.
  * [3] E.W. Barnes.  The theory of the double gamma function.  Phil. Trans. Royal Soc. London (A), 196:265–387., 1901.
  * [4] E.W. Barnes.  On the theory of the multiple gamma functions.  Trans. Cambridge Philos. Soc., 19:374–439., 1904.
  * [5] G. Beylkin and L. Monzón.  On approximation of functions by exponential sums.  Applied and Computational Harmonic Analysis, 19(1):17–48, 2005\.  https://doi.org/10.1016/j.acha.2005.01.003.
  * [6] G. Beylkin and L. Monzón.  Approximation by exponential sums revisited.  Applied and Computational Harmonic Analysis, 28(2):131–149, 2010\.  https://doi.org/10.1016/j.acha.2009.08.011.
  * [7] D. Braess and W. Hackbusch.  Approximation of $1/x$ by exponential sums in $[1,\infty)$.  IMA Journal of Numerical Analysis, 25(4):685–697, 2005.  https://doi.org/10.1093/imanum/dri015.
  * [8] C.-P. Chen.  Asymptotic expansions for Barnes $G$-function.  Journal of Number Theory, 135:36–42, 2014.  https://doi.org/10.1016/j.jnt.2013.08.007.
  * [9] C. Ferreira and J. L. López.  An asymptotic expansion of the double gamma function.  Journal of Approximation Theory, 111(2):298–314, 2001.  https://doi.org/10.1006/jath.2001.3578.
  * [10] J. P. Keating and N. C. Snaith.  Random matrix theory and $\zeta(1/2+it)$.  Communications in Mathematical Physics, 214:57–89, 2000.  https://doi.org/10.1007/s002200000261.
  * [11] A. Kuznetsov.  On extrema of stable processes.  The Annals of Probability, 39(3):1027 – 1060, 2011.
  * [12] C. Lanczos.  A precision approximation of the gamma function.  J. SIAM Numer. Anal. B, 1:86–94, 1964.  https://doi.org/10.1137/0701008.
  * [13] L. Lewin.  Polylogarithms and associated functions.  Elsevier North Holland, New York, 1981.
  * [14] L.C. Maximon.  The dilogarithm function for complex argument.  Proc. R. Soc. Lond. A, 459:2807–2819, 2003.  http://doi.org/10.1098/rspa.2003.1156.
  * [15] G. Nemes.  Error bounds and exponential improvement for the asymptotic expansion of the Barnes $G$-function.  Proceedings of the Royal Society A, 470(2172), 2014.  https://doi.org/10.1098/rspa.2014.0534.
  * [16] A. Nikeghbali and M. Yor.  The Barnes $G$ function and its relations with sums and products of generalized Gamma convolution variables.  Electronic Communications in Probability, 14:396 – 411, 2009.  https://doi.org/10.1214/ECP.v14-1488.
  * [17] G. R. Pugh.  An analysis of the Lanczos gamma approximation.  PhD thesis, The University of British Columbia, 2004.  https://dx.doi.org/10.14288/1.0080001.
  * [18] J. R. Rice.  Chebyshev approximation by exponentials.  Journal of the Society for Industrial and Applied Mathematics, 10(1):149–161, 1962.  http://www.jstor.org/stable/2098804.
  * [19] P. Sarnak.  Determinants of Laplacians.  Communications in Mathematical Physics, 110:113–120, 1987.  https://doi.org/10.1007/BF01209019.
  * [20] T. Shintani.  A Proof of the Classical Kronecker Limit Formula.  Tokyo Journal of Mathematics, 3(2):191 – 199, 1980.  https://doi.org/10.3836/tjm/1270472992.
  * [21] I. Vardi.  Determinants of Laplacians and multiple gamma functions.  SIAM Journal on Mathematical Analysis, 19(2):493–507, 1988.  https://doi.org/10.1137/0519035.
  * [22] M.F. Vignéras.  L’équation fonctionalie de la fonction zeta de Selberg de groupe modulaire $PSL(2,{\mathbb{Z}})$.  Asterisque, 61:235–249, 1979.
  * [23] W. Jedidi, T. Simon and M. Wang.  Density solutions to a class of integro-differential equations.  Journal of Mathematical Analysis and Applications, 458(1):134–152, 2018.  https://doi.org/10.1016/j.jmaa.2017.08.043.
  * [24] Z. Xu and W. Wang.  More asymptotic expansions for the Barnes $G$-function.  Journal of Number Theory, 174:505–517, 2017.  https://doi.org/10.1016/j.jnt.2016.11.007.

##  Appendix A Computing the dilogarithm function for $|z|\leq 1$

We will follow closely the discussion in Section 4 in [14]. An efficient
method for computing ${\textrm{Li}}_{2}(z)$ for $|z|\leq 1$, which avoids the
issue of slow convergence of the series (18) when $|z|$ is close to one, can
be obtained via the so-called Debye function

$D(w):=\int_{0}^{w}\frac{u{\textnormal{d}}u}{e^{u}-1}.$ (A.32)

It is easy to see (by substitution $t=1-e^{-u}$) that

${\textrm{Li}}_{2}(z)=D(-\ln(1-z)).$ (A.33)

The Debye function has a Taylor series expansion

$D(w)=w-\frac{w^{2}}{4}+\sum\limits_{n\geq 1}B_{2n}\frac{w^{2n+1}}{(2n+1)!},$
(A.34)

where $B_{n}$ are Bernoulli numbers. From the asymptotics of Bernoulli numbers

$|B_{2n}|=2(2n)!(2\pi)^{-2n}(1+o(1)),\;\;\;n\to+\infty$

it is easy to deduce that the above Taylor series converges for $|w|<2\pi$
(this fact also follows by noting that the integrand in (A.32) is analytic in
the disk $|u|<2\pi$, which implies that Debye function $D(w)$ must be analytic
in $|w|<2\pi$).

Now, our algorithm for computing ${\textrm{Li}}_{2}(z)$ for $|z|\leq 1$
consists of the following two steps:

  * (i)

If $z\in{\Omega_{1}}:=\\{z\in{\mathbb{C}}\;:\;|z|\leq
1,\;\textnormal{Re}(z)\leq 1/2\\}$ then we use identity (A.33) and compute
$D(-\ln(1-z))$ via series expansion (A.34). One can check that for
$z\in{\Omega}_{1}$ we have $|\ln(1-z)|\leq\pi/3$, thus the terms in the series
(A.34) decay at the rate $O(6^{-2n})$.

  * (ii)

If $z\in{\Omega_{2}}:=\\{z\in{\mathbb{C}}\;:\;|z|\leq
1,\;\textnormal{Re}(z)>1/2\\}$ then we first use identity (3.3) in [14]

${\textrm{Li}}_{2}(z)=-{\textrm{Li}}_{2}(1-z)+\frac{\pi^{2}}{6}-\ln(1-z)\ln(z).$

If $z\in\Omega_{2}$ then $1-z\in\Omega_{1}$, thus ${\textrm{Li}}_{2}(1-z)$ can
be computed by the same approach as described in item (i) above.

##  Appendix B Coefficients $\lambda_{j}$ and $c_{j}$ in the case $p=15$,
$N_{1}=22$ and $N_{2}=8$

Here we present values of $\\{\lambda_{j}\\}_{1\leq j\leq 15}$ and
$\\{c_{j}\\}_{1\leq j\leq 15}$. The complex values come in conjugate pairs,
thus we present only one such value and the other value may be computed by
taking conjugate of the previous number (in other words,
$\lambda_{2n+1}=\overline{\lambda_{2n}}$ and $c_{2n+1}=\overline{c_{2n}}$ for
$4\leq n\leq 7$ ). These lists of numbers can also be downloaded as text files
from the author’s webpage.

$\displaystyle\lambda_{1}$
$\displaystyle=1.015816941860969308,\;\;\;\lambda_{2}=1.053963061918305102,$
$\displaystyle\lambda_{3}$
$\displaystyle=1.116651540074509609,\;\;\;\lambda_{4}=1.207738507792217625,$
$\displaystyle\lambda_{5}$
$\displaystyle=1.332888622825204091,\;\;\;\lambda_{6}=1.719941572880692604,$
$\displaystyle\lambda_{7}$ $\displaystyle=2.930503690937967271,$
$\displaystyle\lambda_{8}$
$\displaystyle=2.231464874614817990-{\textnormal{i}}\times
0.280912039207008020,$ $\displaystyle\lambda_{10}$
$\displaystyle=2.639898812086004465-{\textnormal{i}}\times
0.873853916915943961,$ $\displaystyle\lambda_{12}$
$\displaystyle=2.941124258312725471-{\textnormal{i}}\times
1.605727317761697042,$ $\displaystyle\lambda_{14}$
$\displaystyle=3.229198135526167105+{\textnormal{i}}\times
2.596457178929701727,$ $\displaystyle c_{1}$
$\displaystyle=-3.361986110456561101\times
10^{-5},\;\;\;c_{2}=-1.894144561517152089\times 10^{-4},$ $\displaystyle
c_{3}$ $\displaystyle=-5.010483210821698243\times
10^{-4},\;\;\;c_{4}=-8.578556468220969250\times 10^{-4},$ $\displaystyle
c_{5}$ $\displaystyle=-8.943696088058549902\times
10^{-4},\;\;\;c_{6}=1.854241163038972664\times 10^{-3},$ $\displaystyle c_{7}$
$\displaystyle=-1.918606889602829249\times 10^{-5},$ $\displaystyle c_{8}$
$\displaystyle=-3.849191533344471619\times 10^{-4}+{\textnormal{i}}\times
2.988868248105834482\times 10^{-4},$ $\displaystyle c_{10}$
$\displaystyle=1.121264751590328248\times 10^{-5}-{\textnormal{i}}\times
4.979727219667585924\times 10^{-6},$ $\displaystyle c_{12}$
$\displaystyle=-1.113878636296735895\times 10^{-7}-{\textnormal{i}}\times
9.472403853117676266\times 10^{-8},$ $\displaystyle c_{14}$
$\displaystyle=-1.508505417972961883\times 10^{-10}-{\textnormal{i}}\times
3.899201018438800852\times 10^{-10}.$