# Ophthalmic Biomarker Detection Using Ensembled Vision Transformers - Winning
Solution to IEEE SPS VIP Cup 2023

###### Abstract

This report outlines our approach in the IEEE SPS VIP Cup 2023: Ophthalmic
Biomarker Detection competition. Our primary objective in this competition was
to identify biomarkers from Optical Coherence Tomography (OCT) images obtained
from a diverse range of patients. Using robust augmentations and 5-fold cross-
validation, we trained two vision transformer-based models: MaxViT and EVA-02,
and ensembled them at inference time. We find MaxViT’s use of convolution
layers followed by strided attention to be better suited for the detection of
local features while EVA-02’s use of normal attention mechanism and knowledge
distillation is better for detecting global features. Ours was the best-
performing solution in the competition, achieving a patient-wise F1 score of
0.814 in the first phase and 0.8527 in the second and final phase of VIP Cup
2023, scoring 3.8% higher than the next-best solution.

##  1 Introduction

Accurate ophthalmic biomarker detection using Optical Coherence Tomography
(OCT) images has received tremendous attention in contemporary research in
ophthalmology, with significant implications in the diagnosis and treatment of
eye conditions. In this context, the IEEE SPS VIP Cup 2023 presented a unique
opportunity to explore innovative approaches in the field of ophthalmic
biomarker detection. This paper outlines our approach to address this issue,
focusing on utilizing vision transformer-based models to analyze OCT images.

On this backdrop, to solve this ophthalmic biomarker detection task, we
employed two distinct models, the Multi-Axis Vision Transformer (MaxViT) [1]
and EVA-02 [2], carefully selected following a systematic exploration, for
their respective proficiency in identifying local and global features within
the images. The intricate exploration and subsequent integration of these
models were pivotal in formulating a solution that not only fulfilled the
competition’s criteria but also advanced the state-of-the-art of ophthalmic
biomarker detection in both understanding and methodology.

##  2 Methodology

###  2.1 Dataset

The competition provided a rich dataset, OLIVES [3], encompassing 9408 labeled
image-biomarker pairs from 96 patients and an additional 78185 unlabeled OCT
images, each accompanied by clinical labels. The competition was divided into
two phases as follows. In Phase 1 (2), the test dataset consisted of 3871
(250) images from 40 different (167 new) patients.

Each OCT scan segment was labeled to denote the presence or absence of 6
biomarkers, namely Intraretinal Hyperreflective Foci (IRHRF), Partially
Attached Vitreous Face (PAVF), Fully Attached Vitreous Face (FAVF),
Intraretinal Fluid (IRF), Diffuse Retinal Thickening or Diabetic Macular Edema
(DRT/DME) and Vitreous Debris (VD). Depending on the spatial extent, IRHRF and
IRF can be loosely grouped as local features, meaning they could be detected
by looking at just a subsection of the image. On the other hand, PAVF, FAVF,
and VD are global features with DRT/DME falling somewhat in between. We
elucidate the rationale behind the dataset partitioning in relation to model
architecture in section 3.3.

###  2.2 Models Considered

We considered multiple variants of ResNet [4] models and Inception [5] models
(collectively referred to as Convolution-based Models henceforth).

Inspired by [6], we added Convolutional Block Attention Modules (CBAM) [7] to
InceptionResnetV2 (referred to as IRV2_CBAM for brevity). We added three such
CBAMs after the Stem, Reduction A, and Reduction B modules of
InceptionResnetV2. The improved performance of IRV2_CBAM (to be presented in
Section 3) inspired us to move to vision transformer models, including ViT
[8], MaxViT [1], and EVA-02 [2].

Our early tests indicated an important role for image dimensions when
detecting biomarkers. This observation was corroborated through a consultation
with an ophthalmologist, wherein the discussion came up that downsizing images
to a resolution of 224x224 pixels might have made it harder to identify these
biomarkers. As such, we focused on models pre-trained on larger images. ViT
[8], MaxViT [1] and EVA-02 [2] support image resolutions of $384\times 384$,
$512\times 512$ and $448\times 448$ respectively. Notably, we could only use
the base version of these models due to computational constraints.

###  2.3 Hyperparameters

We used AdamW[9] optimizer with default initialization and set the initial
learning rate to $3\times 10^{-5}$. We used the Exponential Learning Rate
Scheduler, with a weight decay of 0.9. For convolution-based models, we used
128 as the batch size and trained models for 35 epochs, with early stopping
based on the best cross-validation F1 score. For transformer-based models, we
used the largest possible batch size supported by our hardware, which was 1
for MaxViT and 2 for both EVA-02 and ViT. To account for the small batch size,
we set the gradient accumulation steps to 8. We trained all vision transformer
models for two epochs. We found all ViT models to be prone to overfitting the
training data after 2 epochs.

###  2.4 Data Augmentation

In Phase 1, we used random greyscale transformation with $p=0.2$, color jitter
with $p=0.8$, random resized crop with $scale=(0.7,1)$, random horizontal
flip, and finally, normalization with a mean of 0.1706 and a standard
deviation of 0.2112. We found 0.7 to be the optimal scale for random resized
crop while keeping other augmentations constant. Other augmentation parameters
were not optimally tuned.

For Phase 2, we add a random perspective shift augmentation with $distortion\
scale=0.23$, $p=1$, and $fill=255$ to make the training data similar to the
Phase 2 evaluation dataset. In both phases, we did not augment the test data
beyond resizing and normalization.

###  2.5 5-fold Cross Validation

For both phases, we performed 5-fold cross-validation where we partitioned the
data into 5 folds with 80% in the train set and 20% on the validation set. On
these 5 different folds, we trained our models and ran inference on the test
set after every epoch, and combined the confidence scores to obtain the final
binary decision for each biomarker.

Fig. 1: Our final training pipeline. For evaluation, all augmentations except
Resize were removed.

###  2.6 Ensembling MaxViT and EVA-02

The complementary strengths of MaxViT and EVA-02 naturally imply that
ensembling their outputs has the potential to improve upon their individual
performance across all biomarkers. One straightforward way to implement this
is by using MaxViT to detect local biomarkers and using EVA-02 for global
biomarkers. In this scheme, MaxViT’s predictions for global biomarkers are
entirely ignored (as well as EVA-02’s predictions for local ones). We also
apply a finer-grained ensembling scheme, where we average both model’s output
probabilities. Fig. 1 presents a schematic overview of our overall pipeline.
We will refer to this (finer-grained) ensemble as MaxViT-EVA-02.

###  2.7 Evaluation Metrics

In the domain of medical imaging where severe class imbalance is the norm, the
F1 score often is the metric of choice instead of accuracy. For Phase 1 of the
competition, to test the generalization ability of solutions, the F1 score was
calculated over all the images in the test set. For Phase 2, to measure
personalization: how well a model performs on individual patients, patient-
wise F1 scores were calculated over images from the same patient and these
scores were averages over all patients in the test dataset. More details can
be found on the competition website.

###  2.8 Hardware Specification and Environment Setup

For convolution-based models implemented in Tensorflow[10], we used Kaggle TPU
VM v3-8 instances paired with 330GB RAM. Due to the limited support of state-
of-the-art models on TPU, we mainly used this setup for pilot experiments. For
transformer-based models (implemented in PyTorch 2.0.1[11] and ‘timm’ [12]
library with the weights hosted on Hugging Face), we used Kaggle Nvidia P100
GPU instances with 16GB VRAM, 13GB RAM, and 19GB disk space. We used scikit-
learn [13] libraries for other auxiliary needs. The runtime of our complete
MaxViT pipeline, including training, validation, and inference, was
approximately 11 hours, while that of our EVA-02 pipeline was approximately 7
hours.

##  3 RESULTS AND DISCUSSIONS

###  3.1 Model Selection Results

To establish a baseline, we trained multiple variants of ResNet [4] models and
Inception [14] models. We find that model size or model performance on
ImageNet dataset [15] are not reliable indicators of its suitability for the
task at hand (Table 1). InceptionResnetV2[5] (55.84 M parameters) proved to be
the most effective model with an F1 score of 0.686 and the much smaller
InceptionV3 (23.83 M parameters) model performed comparably with an F1 score
of 0.682 (Table 1).

Model | Param(M) | ImageNet | Test F1
---|---|---|---
ConvNextBase | 88.59 | 87.13 | 0.612
Resnet50 | 25.57 | 75.30 | 0.634
Resnet152 | 66.84 | 78.57 | 0.649
Resnet101 | 44.57 | 78.25 | 0.657
EfficientNetV2L | 118.52 | 86.80 | 0.662
InceptionV3 | 23.83 | 78.95 | 0.682
InceptionResnetV2 | 55.84 | 80.46 | 0.686
Table 1: Comparison of Convolution-based Models. We report the number of model
parameters, Top1 Accuracy on the ImageNet [15] dataset collected from
PapersWithCode, and F1 score on the Phase 1 Test dataset. All models were
evaluated using 5-fold cross-validation.

###  3.2 5-fold Cross Validation Results

5-fold cross-validation boosts out Phase-1 test scores substantially. Initial
experiments revealed that our best-performing convolution-based model,
InceptionResnetV2 consistently scored 0.66 when trained on random 80% splits
of the train set. However, using cross-validation, InceptionResnetV2
consistently scored around 0.68. As such, we used cross-validation in all
further experiments and in the final submission as well. Individually, MaxViT
and EVA-02 models scored 0.68 while with cross-validation they scored 0.71.

###  3.3 Classification of Biomarkers according to Spatial Extent

Upon reviewing the images, we noticed that biomarkers B1 (IRHRF), B4 (IRF),
and B5 (DRT/DME) were local, meaning they could be detected by looking at just
a subsection of the image. This observation was confirmed by an
ophthalmologist, who also mentioned that B5 is somewhat in between local and
global.

###  3.4 Analysis Of Adding CBAM

Adding CBAM [7] to InceptionResnetV2 substantially boosted the F1 score from
0.686 to 0.696 (Table 2) for a negligible increase in the network complexity
(i.e., parameter count increased by only 0̃.37%; not reported in the table).
Notably, this boost in performance actually inspired us to move to vision
transformer models.

To understand the reason for the improved F1 scores, we calculated the F1
score across biomarker types individually and discovered that CBAM improved
the performance on certain biomarkers substantially while showing marginal
improvement in others. It even registered a deterioration, albeit only
slightly, in one case. Therefore, we hypothesize that the attention module
improved the detection of local biomarkers.

Biomarker | Type | IRV2 | IRV2_CBAM | VIT_BASE
---|---|---|---|---
IRHRF | L | 0.709 | 0.746 (+) | 0.773 (+)
PAVF | G | 0.610 | 0.609 | 0.662 (+)
FAVF | G | 0.837 | 0.841 | 0.869 (+)
IRF | L | 0.557 | 0.599 (+) | 0.552
DRT/DME | L/G | 0.599 | 0.628 (+) | 0.594
VD | G | 0.753 | 0.759 | 0.755
Overall |  | 0.686 | 0.696 | 0.701
Table 2: Comparison of InceptionResnetV2 with (IRV2_CBAM) and without (IRV2)
CBAM. L (G) in the type column refers to Local (Global). For individual
biomarker types, a plus sign in the bracket beside a score indicates
significant improvement against the score of the network to its immediate left
column. All models were evaluated using 5-fold cross-validation.

###  3.5 Comparison of Convolution-based and Transformer-based Models

Although adding an attention mechanism in the form of CBAM to InceptionResnet
specifically improves the performance on local biomarkers, we find no such
correlation when comparing convolution-based models and the purely attention-
based ViT [8] architectures. This suggests the need for explicit convolution
in addition to attention for optimal biomarker detection.

###  3.6 Effectiveness of Convolution and Attention

MaxViT[1] is a vision transformer model composed of multiple MaxViT blocks
where each block performs convolution, strided/block attention, and
dilated/grid attention. The addition of explicit convolution makes MaxViT
ideal for biomarker detection. We achieved an F1 score of 0.718 (Table 3)
using the base variant of the MaxViT model, which is a substantial improvement
over IRV2_CBAM and ViT_BASE. However, MaxViT does not utilize true attention
across all image tokens, which motivated us to test EVA-02 [2], a plain Vision
Transformer model that improves upon the standard ViT [8] by using a 1B
parameter EVA-CLIP model as its teacher. The parameter counts of MaxViT and
EVA-02 are 119.88M and 87.12M respectively. Comparing MaxViT and EVA-02 across
the 6 biomarkers, we see that EVA-02 performs noticeably better on global
biomarkers despite being smaller between the two. we hypothesize that MaxViT’s
sparse attention improves the detection of local biomarkers while EVA-02’s
true attention excels at detecting global features.

###  3.7 Ensembling Results

While our simple ensembling does boost the test set F1 score to 0.720 (not
shown in the table for brevity), the finer-grained ensembling scheme yields an
even greater performance with an improved F1 score of 0.724.

Biomarker | Type | MaxViT | EVA-02 | Ensemble
---|---|---|---|---
IRHRF | L | 0.774 | 0.731 | 0.779
PAVF | G | 0.677 | 0.701 | 0.688
FAVF | G | 0.868 | 0.874 | 0.879
IRF | L | 0.611 | 0.575 | 0.600
DRT/DME | L/G | 0.615 | 0.593 | 0.618
VD | L | 0.764 | 0.779 | 0.782
Overall | Overall | 0.718 | 0.709 | 0.724
Table 3: Comparison of MaxViT, EVA-02, and their Ensemble for Various
Biomarkers. The models have been ensembled by averaging their output
probabilities.(L: Local, G: Global)

###  3.8 Patient-wise F1 score and Leaderboard Position

Our MaxViT-EVA-02 ensemble pipeline achieved a patient-wise F1 score of 0.814
averaging over 40 patients and 3781 images in the first phase and 0.8527 in
the second phase (167 patients and 250 images). Our second phase F1 score is
3.8% higher than the next best solution (0.8215) in the competition as per the
competition website.

###  3.9 Leveraging Unlabelled Training Data

Several strategies were explored to enhance the model’s performance on the
task using unlabelled data and the accompanying clinical labels. Initially,
contrastive learning was applied following [16]; however, it did not yield any
improvement, maintaining a score of 0.686. Although looked promising, we were
not able to explore contrastive learning using ViT models due to computational
constraints. Subsequently, pretraining the InceptionResnetV2 model was
attempted, but this approach also failed to enhance the model’s performance,
with the score remaining stagnant at 0.686. To exploit the availability of two
additional clinical labels besides the six target biomarkers, the model was
trained to predict all eight labels, anticipating a more refined gradient
descent. Unfortunately, this modification did not lead to any improvement,
with the score persisting at 0.686. Lastly, the implementation of pseudo-
labeling was explored, but it significantly deteriorated the model’s
performance, plummeting the score to 0.519.

###  3.10 Analysis of Outlying Patient-wise F1 Scores

In the analysis of cases where the model exhibited a low F1 score in detecting
biomarkers from OCT scans, several patterns were observed. Patient 01-002 at
week 40 and patient 02-044 at week 0 presented with severe spots, resulting in
F1 scores of 0.64 and 0.55, respectively. Moderate spots were identified as
the likely cause for the low F1 scores of 0.6 in patients 01-007 at week 100
and 01-049 at week 0. Additionally, patient 01-043 at week 100 exhibited a
severe artifact, leading to the lowest F1 score of 0.37. Moderate artifacts
were also noted in patients 01-049 and 02-044 at week 100, with F1 scores of
0.6 and 0.52, respectively. However, the likely cause for the low F1 scores
observed in patients 01-019, 01-036, and 01-054 at week 100 (F1 scores of
0.51, 0.62, and 0.48) are not immediately evident to non-medical
professionals. We leave a more thorough analysis and subsequent pipeline
adjustments as future work.

##  4 CONCLUSION

In this work, we presented the methodology of our final submission to the
second phase of IEEE SPS VIP Cup 2023: Ophthalmic Biomarker Detection. We also
presented the underlying motivation for pipeline design decisions. We find
that Vision Transformer (ViT) models have begun to consistently outperform
their Convolutional Neural Network (CNN) counterparts. Furthermore, we find
that k-fold cross-validation and model ensembling continue to be effective
means to utilize the entire dataset and to improve the generalization of
predictions.

##  5 ACKNOWLEDGEMENT

We would like to extend our sincere gratitude to Dr. S.M. Rezwan Hussain, a
distinguished ophthalmologist at the Eye Department, Combined Military
Hospital(CMH), Dhaka. Bangladesh, for his invaluable insights and expertise
regarding biomarker classification according to their spatial extent.

## References

  * [1] Zhengzhong Tu, Hossein Talebi, Han Zhang, Feng Yang, Peyman Milanfar, Alan Bovik, and Yinxiao Li,  “Maxvit: Multi-axis vision transformer,”  ECCV, 2022.
  * [2] Yuxin Fang, Quan Sun, Xinggang Wang, Tiejun Huang, Xinlong Wang, and Yue Cao,  “Eva-02: A visual representation for neon genesis,”  arXiv preprint arXiv:2303.11331, 2023.
  * [3] Mohit Prabhushankar, Kiran Kokilepersaud, Yash-yee Logan, Stephanie Trejo Corona, Ghassan AlRegib, and Charles Wykoff,  “Olives dataset: Ophthalmic labels for investigating visual eye semantics,”  2022\.
  * [4] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun,  “Deep residual learning for image recognition,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 770–778.
  * [5] Christian Szegedy, Sergey Ioffe, Vincent Vanhoucke, and Alexander Alemi,  “Inception-v4, inception-resnet and the impact of residual connections on learning,”  in Proceedings of the AAAI conference on artificial intelligence, 2017, vol. 31.
  * [6] Zhuang Ai, Xuan Huang, Jing Feng, Hui Wang, Yong Tao, Fanxin Zeng, and Yaping Lu,  “Fn-oct: Disease detection algorithm for retinal optical coherence tomography based on a fusion network,”  Frontiers in Neuroinformatics, vol. 16, 2022.
  * [7] Sanghyun Woo, Jongchan Park, Joon-Young Lee, and In So Kweon,  “Cbam: Convolutional block attention module,”  in Proceedings of the European conference on computer vision (ECCV), 2018, pp. 3–19.
  * [8] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, et al.,  “An image is worth 16x16 words: Transformers for image recognition at scale,”  arXiv preprint arXiv:2010.11929, 2020.
  * [9] Ilya Loshchilov and Frank Hutter,  “Decoupled weight decay regularization,”  arXiv preprint arXiv:1711.05101, 2017.
  * [10] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng,  “TensorFlow: Large-scale machine learning on heterogeneous systems,” 2015,  Software available from tensorflow.org.
  * [11] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala,  “Pytorch: An imperative style, high-performance deep learning library,”  in Advances in Neural Information Processing Systems 32, pp. 8024–8035. Curran Associates, Inc., 2019.
  * [12] Ross Wightman,  “Pytorch image models,” https://github.com/rwightman/pytorch-image-models, 2019.
  * [13] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al.,  “Scikit-learn: Machine learning in python,”  Journal of machine learning research, vol. 12, no. Oct, pp. 2825–2830, 2011.
  * [14] Christian Szegedy, Wei Liu, Yangqing Jia, Pierre Sermanet, Scott Reed, Dragomir Anguelov, Dumitru Erhan, Vincent Vanhoucke, and Andrew Rabinovich,  “Going deeper with convolutions,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2015, pp. 1–9.
  * [15] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei,  “Imagenet: A large-scale hierarchical image database,”  in 2009 IEEE conference on computer vision and pattern recognition. Ieee, 2009, pp. 248–255.
  * [16] Kiran Kokilepersaud, Stephanie Trejo Corona, Mohit Prabhushankar, Ghassan AlRegib, and Charles Wykoff,  “Clinically labeled contrastive learning for oct biomarker classification,”  IEEE Journal of Biomedical and Health Informatics, 2023.


# On the non-vanishing of modular $L$-values and Fourier coefficients of cusp
forms

Jun hwi Min<EMAIL_ADDRESS>Ulsan National Institute of Science and
Technology, Ulsan, Korea

###### Abstract.

We prove a non-vanishing result of modular $L$-values with quadratic twists,
where the quadratic discriminants are in a short interval. Using this fact and
Waldspurger’s theorem, we improve the results of Balog-Ono[The chebotarev
density theorem in short intervals and some questions of Serre, Journal of
number theory. 91(2):356-371(2001)] on the non-vanishing of Fourier
coefficients of half-integral weight eigenform.

###### Key words and phrases:

Modular $L$-values, Fourier coefficients, Non-vanishing

##  1\. Introduction and statement of results

Let $f$ be a Hecke eigenform of level $N$ and weight $2k$ where $k$ is a
positive integer, and let $\chi_{d}$ be a primitive quadratic character
associated to the quadratic field $\mathbb{Q}(\sqrt{d})$. The study of the
non-vanishing of $L(k,f,\chi_{d})$ arose from the study of the rank of modular
elliptic curves. The celebrated theorem of Kolyvagin asserted that if $L(s,E)$
has a simple zero at $s=1$ and there exists a quadratic discriminant $d<0$
coprime to $4N$ such that $L(1,E,\chi_{d})\neq 0$, then the
rank$E(\mathbb{Q})=1$(See the introduction of [13]). Murty-Murty[9] showed
that one can get rid of the non-vanishing condition, by showing that there
exists an infinite family of quadratic discriminants with $L(1,E,\chi_{d})\neq
0$. Even further, the quantitative estimates on the number of quadratic
discriminants with $L(1,E,\chi_{d})\neq 0$ has gained huge interests. Let

$\displaystyle N_{f}(X):=\\#\\{d:d\text{ is a square-free quadratic
discriminant},|d|\leq X,L(k,f,\chi_{d})\neq 0\\}.$

The well-known conjecture on the non-vanishing of $L(k,f,\chi_{d})$ asserts
that $N_{f}(X)\gg X$. The best known lower bound is that of Ono-Skinner[10];
when $f$ is a even weight newform of trivial nebentype,

$\displaystyle N_{f}(X)\gg X/\log X.$

We are interested in a short interval analogue of $N_{f}(X)$, where $f(z)\in
S_{2k}(N,\psi)$ is a normalized Hecke eigenform. Here, the set of quadratic
discriminants are defined modulo $4N$, Namely

$\displaystyle\mathcal{D}:=\\{0<(-1)^{k}d:d\equiv v^{2}\mod 4N\text{ for some
}(v,4N)=1\\}.$

We denote $\gamma(4N):=\\#\mathcal{D}$. The counting function of the number of
nonzero $L$-values in short interval is defined by

$\displaystyle N_{f}(X,h):=\\#\\{d\in\mathcal{D}:d\text{ square-
free},X\leq|d|\leq X+h,L(k,f,\chi_{d})\neq 0\\}.$

To bound $N_{f}(X,h)$, it is necessary to estimate the first moment of
$L$-values in a short interval. That is,

$\displaystyle S_{f}(X,h):$
$\displaystyle=\sideset{}{{}^{\prime}}{\sum}_{\begin{subarray}{c}d\in\mathcal{D}\\\
\ X\leq|d|\leq X+h\end{subarray}}L(k,f,\chi_{d}).$

Here, the superscript ′ denotes that the sum is taken over the square-free
numbers $d$.

######  Theorem 1.1.

For $X^{3/4+\epsilon}\leq h\leq X$, we have

$\displaystyle S_{f}(X,h)=C_{N}L_{f}(k)h+O_{f,\epsilon}(hX^{-\epsilon}),$

where the constants $C_{N}$ and $L_{f}(k)$ are given by

$\displaystyle C_{N}=\frac{3\gamma(4N)}{\pi^{2}N}\prod_{p|4N}(1-p^{-2})^{-1}$

and

(1.1) $\displaystyle
L_{f}(k)=\sum_{\begin{subarray}{c}n=rj^{2}\end{subarray}}\frac{a_{n}}{n^{k}}\prod_{\begin{subarray}{c}p|n\\\
(p,4N)=1\end{subarray}}(1+p^{-1})^{-1}\neq 0,$

respectively. In the R.H.S of (1.1), the indices $r$ and $j$ are positive
integers with $r|(4N)^{\infty}$ and $(j,4N)=1$, respectively. Note $L_{f}(s)$
is the $L$-function $L_{f,1}(s)$ given in [8, p.385].

To obtain Theorem 1.1, we need to make use of the second moment of $L$-values,
given in [12] and [5].

######  Theorem 1.2.

[12, 5] If $f$ is an eigenform, then

$\displaystyle\sideset{}{{}^{\prime}}{\sum}_{\begin{subarray}{c}|d|\leq
X\end{subarray}}|L(k,f,\chi_{d})|^{2}\ll_{\epsilon}X^{1+\epsilon}.$

The authors in [12] and [5] proved Theorem 1.2 for newforms, but the
discussion goes through for eigenforms as well.

By Cauchy’s inequality, we have an immediate corollary.

######  Corollary 1.3.

For $X^{3/4+\epsilon}\leq h\leq X$, we have

$N_{f}(X,h)\gg\frac{h^{2}}{X^{1+\epsilon}}.$

The implied constant only depends on $f,\epsilon$.

As far as we know, the only present result on non-vanishing of
$L(k,f,\chi_{d})$ in short intervals other than Corollary 1.3 is that of
Balog-Ono[1]. Their result is as follows.

######  Theorem 1.4.

[1] Let $f\in S_{2k}^{new}(\Gamma_{0}(N),\psi_{triv})$. If $f(z)$ is not a
linear combination of weight-$3/2$ theta functions, then there exists a
positive integer $k_{f}$ such that for $X^{1-1/k_{f}+\epsilon}\leq h\leq X$,

(1.2) $\\#\\{d:d\text{ square-free},X\leq|d|\leq X+h,L(k,f,\chi_{d})\neq
0\\}\gg\frac{h}{\log X}.$

Note that in (1.2), $d$ runs through not only $\mathcal{D}$ but also all
square-free quadratic discriminants. Unlike Theorem 1.4, Corollary 1.3 gives a
non-vanishing result for the eigenforms of non-trivial nebentypus. In case of
trivial nebentype, although Corollary 1.3 gives a weaker result when $h\geq
X^{1-1/k_{f}+\epsilon}$, it is worth noting that it gives further information
on $N_{f}(X,h)$ when $X^{3/4+\epsilon}\leq h<X^{1-1/k_{f}+\epsilon}$. The
advantage of having a non-vanishing result in a shorter inteval is more
visible in view of the non-vanishing of Fourier coefficients.

Theorem 1.4 is a consequence of the following corollary and Waldspurger’s
theorem.

######  Theorem 1.5.

[1, Corollary 4] Let $g\in S_{k+1/2}(M,\psi)$, i.e., a cuspform of level $M$
and weight $k+1/2$ with the nebentype $\psi$. If $g(z)$ is not a linear
combination of weight-$3/2$ theta functions, there exists a positive integer
$k_{f}$ such that for $X^{1-1/k_{f}+\epsilon}\leq h\leq X$,

$\\#\\{X\leq n\leq X+h:a_{g}(n)\neq 0\\}\gg\frac{h}{\log X}.$

######  Theorem 1.6.

[11] Let $N$ be even, $f\in S_{2k}(N,\psi^{2})$ and $g\in S_{k+1/2}(2N,\psi)$
where $f$ is the Shimura correspondent of $g$. Then

$a_{g}(d)^{2}=\kappa_{f}L(k,f,\psi_{0}^{-1}\chi_{d})\psi(d)d^{k},$

where $\psi_{0}(n)=\psi(n)(\frac{-1}{n})^{k}$.

Even further, Theorem 1.5 answers to a stronger form of Serre’s questions.
Serre initiated the study on the sizes of gaps between non-zero Fourier
coefficients. He defined the following gap function.

$\displaystyle i_{f}(n)=\begin{cases}\max\\{i:a_{f}(n+j)=0\text{ for all
}0\leq j\leq i\\}&\text{if }a_{f}(n)=0,\\\ 0&\text{otherwise}.\end{cases}$

Serre’s question asks us to estimate this gap function with regards to the
cusp forms of integral or half-integral weights. Combining Theorem 1.6 and
Corollary 1.3, our second corollary immediately follows.

######  Corollary 1.7.

Assume the conditions in Theorem 1.6 with $\psi=\psi_{triv}$. Then for
$X^{3/4+\epsilon}\leq h\leq X$,

$\displaystyle\\#\\{X\leq n\leq X+h:a_{g}(n)\neq
0\\}\gg_{f,\epsilon}\frac{h^{2}}{X^{1+\epsilon}}.$

In particular, $i_{g}(n)\ll n^{3/4+\epsilon}$.

Corollary 1.7 along with Theorem 1.5 are the only results on $i_{g}(n)$ where
$g$ is of half-integral weight. On the other hand, many mathematicians have
studied $i_{f}(n)$ where $f$ is of integral weight. Unfortunately, most of the
approaches in the study of integral weight case are not available in half-
integral weight case. Let us briefly review the works on $i_{f}(n)$.

Rankin-Selberg estimates, the multiplicativity of Hecke operators, or the
arithmetic of Galois representations have been useful tools to study
$i_{f}(n)$ for the eigenforms of integral weight. For example, we have a
classical result on $L(f\otimes\overline{f})$ such that there exists an
integer $c_{f}$ for which

$\displaystyle\sum_{n\leq
X}|a_{f}(n)|^{2}n^{1-\frac{k}{2}}=c_{f}X+O(X^{\frac{3}{5}}).$

It immediately follows that $i_{f}(n)\ll n^{3/5}$. Another approach appeals to
the theory of $\mathfrak{B}$-free numbers, especially when $f$ is not of CM
type. Let $\mathfrak{B}=\\{b_{i}\\}$ be a set of integers such that

$\displaystyle\sum_{b\in\mathfrak{B}}\frac{1}{b}<\infty\text{ and
}(b_{i},b_{j})=1\text{ whenever }i\neq j.$

We say that a natural number $n$ is $\mathfrak{B}$-free if it is not divisible
by any of the elements of $\mathfrak{B}$. Specifically, one can define
$\mathfrak{B}$ as follows.

$\displaystyle\mathfrak{B}=\\{p\text{ prime }:a_{f}(p)=0\\}\cup\\{p|N\text{
prime}\\}.$

$\mathfrak{B}$ has a zero density due to the following result of Serre[14,
p.174, Cor.2]. Let $f$ be a newform with integral weight $2k\geq 2$ which is
not of CM type. Then

(1.3) $\displaystyle\\#\\{p\leq X\text{ prime
}:a_{f}(p)=0\\}\ll_{f,\epsilon}\frac{X}{(\log X)^{3/2-\epsilon}}.$

In view of (1.3) and the multiplicativity of Hecke eigenvalues, $a_{f}(n)\neq
0$ if $n$ is square-free and $\mathfrak{B}$-free. Thus estimating $i_{f}(n)$
becomes a problem of counting $\mathfrak{B}$-free numbers in the short
intervals. Balog-Ono were the first to take this approach and they deduced
that

$\displaystyle i_{f}(n)\ll n^{17/41+\epsilon}.$

Later, the estimates of $i_{f}(n)$ has been refined several times. The best
bound for $i_{f}(n)$ is due to Kowalski-Robert-Wu[6]. They proved that for any
holomorphic non-CM cuspidal eigenform $f$ on general congruence groups,

$\displaystyle i_{f}(n)\ll n^{7/17+\epsilon}.$

If $f$ is of CM type, there are no similar general results on $i_{f}(n)$. The
major difficulty in this case is that the density estimate in (1.3) is valid
only for non-CM type forms. Thus the previous works on CM type forms took
different approaches other than distribution of $\mathfrak{B}$-free numbers.
For example, Das-Ganguly[2] showed that for all nonzero cuspforms $f$ of level
one,

(1.4) $\displaystyle i_{f}(n)\ll n^{1/4}.$

The main ingredient of their work were the congruence relation of Hecke
eigenvalues due to Hatada[3] and the distribution of sum of two squares in
short intervals.

(1.4) can be extended to the eigenforms of higher levels, under some
conditions on $f$. Let $E/\mathbb{Q}$ be an elliptic curve which has a cyclic
rational $4$-isogeny and $f_{E}$ be a newform corresponding to $E$ by the
modularity theorem. Kumar[7] first proved that $f_{E}$ satisfies (1.4), by
showing that there exists a positive integer $m$ such that $a_{f_{E}}(m)\neq
0$ and $m$ is a sum of two squares in intervals $(X,X+cX^{\frac{1}{4}})$. From
this, he deduced that if $f$ is $2$-adically close enough to $f_{E}$, $f$ also
satisfies (1.4).

Let us recall some necessary facts for the future discussion. Let
$e(z)=e^{2\pi iz}$ and

$\displaystyle f(z)=\sum_{n\geq 1}a_{n}e(nz)$

be the Fourier expansion of $f$ at the cusp $\infty$. The associated modular
$L-$function $L(s,f)=\sum_{n\geq 1}a_{n}n^{-s}$ absolutely converges in
$\Re(s)>k+1/2$. It can be analytically continued to an entire function and
satisfies the following functional equation

$\bigg{(}\frac{\sqrt{N}}{2\pi}\bigg{)}^{s}\Gamma(s)L(s,f)=w\bigg{(}\frac{\sqrt{N}}{2\pi}\bigg{)}^{2k-s}\Gamma(2k-s)L(2k-s,f),$

where $w=\pm 1$. In addition, the twisted $L-$function $L(s,f,\chi_{d})$
satisfies the functional equation

$\bigg{(}\frac{d\sqrt{N}}{2\pi}\bigg{)}^{s}\Gamma(s)L(s,f,\chi_{d})=w_{d}\bigg{(}\frac{d\sqrt{N}}{2\pi}\bigg{)}^{2k-s}\Gamma(2k-s)L(2k-s,f,\chi_{d}),$

where $w_{d}=w\chi_{d}(-N)=1$ for all $d\in\mathcal{D}$.

If $f$ is an eigenform, we have the Euler product

(1.5) $\displaystyle
L(s,f)=\prod_{p}\bigg{(}1-\frac{\alpha_{p}}{p^{s}}\bigg{)}^{-1}\bigg{(}1-\frac{\beta_{p}}{p^{s}}\bigg{)}^{-1}$

and it satisfies the Ramanujan-Petersson conjecture (i.e.
$|\alpha_{p}|=|\beta_{p}|=p^{k-1/2}$ for all $p\nmid N$ and
$|\alpha_{p}|,|\beta_{p}|\leq p^{k-1/2}$ otherwise). Given
$\alpha_{p},\beta_{p}$ for all $p$, we define
$n\mapsto\alpha_{n},n\mapsto\beta_{n}$ as totally multiplicative functions on
$\mathbb{N}$.

##  2\. Proof of Theorem 1.1

We modify the methods of [4] and [8]. Necessary changes will be described in
detail.

Set

$\displaystyle V(x)$ $\displaystyle:=\frac{1}{2\pi
i}\int_{(4/5)}\frac{\Gamma(k+s)}{\Gamma(k)}x^{-s}\frac{ds}{s}.$

Here, note that the integral $\int_{(4/5)}$ denotes
$\int_{4/5-i\infty}^{4/5+i\infty}$. By the Mellin transform and the
integration by parts, we have

$\displaystyle
V(x)=\frac{1}{\Gamma(k)}\int_{x}^{\infty}e^{-y}y^{k-1}dy=(1+x+\dots+\frac{x^{k-1}}{(k-1)!})e^{-x}.$

Next, we will define $L(k,f,\chi_{d})$ in terms of the rapidly convergent
sums. Let

$\displaystyle A(Q,\chi_{d})$ $\displaystyle=\frac{1}{2\pi
i}\int_{(4/5)}L(f,\chi_{d},k+s)\frac{\Gamma(k+s)}{\Gamma(k)}\bigg{(}\frac{2\pi}{Q}\bigg{)}^{-s}\frac{ds}{s}.$

We have

$\displaystyle A(Q,\chi_{d})=\sum_{n\geq
1}a_{n}n^{-k}\chi_{d}(n)V\bigg{(}\frac{2\pi n}{Q}\bigg{)}.$

By shifting the contour of integration to $\Re(s)=-4/5$, we obtain

(2.1) $L(k,f,\chi_{d})=A(Q,\chi_{d})+A(d^{2}NQ^{-1},\chi_{d}),$

for any $Q>0$ and square-free $d\in D$. In particular,

(2.2) $L(k,f,\chi_{d})=2A(|d|\sqrt{N},\chi_{d}).$

By Abel’s summation formula,

(2.3) $A(Q,\chi_{d})\ll_{f}Q^{\frac{1}{2}}.$

Combining (2.1) and (2.3) gives

(2.4) $L(k,f,\chi_{d})=A(Q,\chi_{d})+O_{f}(|d|Q^{-\frac{1}{2}})\text{ for all
}Q>0.$

As in [8], we have an upper bound of the fourth moment of $L$-values

$\sideset{}{{}^{\prime}}{\sum}_{d\in D,|d|\leq
X+h}|L(k,f,\chi_{d})|^{4}\ll(X+h)^{2+\epsilon}.$

The first moment we are considering is

$S_{f}(X,h):=\sideset{}{{}^{\prime}}{\sum}_{\begin{subarray}{c}d\in D\\\
X\leq|d|\leq
X+h\end{subarray}}L(k,f,\chi_{d})=2\sideset{}{{}^{\prime}}{\sum}_{\begin{subarray}{c}d\in
D\\\ X\leq|d|\leq X+h\end{subarray}}A(|d|\sqrt{N},\chi_{d}).$

We introduce the Möbius function to relax the square-free condition, so that

$\displaystyle S_{f}(X,h)$
$\displaystyle=2\sum_{\begin{subarray}{c}d\in\mathcal{D}\\\
X/a^{2}\leq|d|\leq(X+h)/a^{2}\end{subarray}}\sum_{a^{2}|d}\mu(a)A(|d|\sqrt{N},\chi_{d})$
$\displaystyle=2\sum_{(a,4N)=1}\sum_{\begin{subarray}{c}d\in\mathcal{D}\\\
X/a^{2}\leq|d|\leq(X+h)/a^{2}\end{subarray}}\mu(a)A(a^{2}|d|\sqrt{N},\chi_{a^{2}d}).$

Here, $A$ is a large number only dependent on $X$. We choose $A$ later. We
split the sum into two parts, say, $S_{f}(X,h)=S+R$, where

$\displaystyle S$ $\displaystyle=2\sum_{\begin{subarray}{c}a\leq A\\\
(a,4N)=1\end{subarray}}\mu(a)\sum_{\begin{subarray}{c}d\in\mathcal{D}\\\
X/a^{2}\leq|d|\leq(X+h)/a^{2}\end{subarray}}A(a^{2}|d|\sqrt{N},\chi_{a^{2}d})$

and

$\displaystyle R$ $\displaystyle=2\sum_{\begin{subarray}{c}a>A\\\
(a,4N)=1\end{subarray}}\mu(a)\sum_{\begin{subarray}{c}d\in\mathcal{D}\\\
X/a^{2}\leq|d|\leq(X+h)/a^{2}\end{subarray}}A(a^{2}|d|\sqrt{N},\chi_{a^{2}d})$
$\displaystyle=2\sum_{(b,4N)=1}\bigg{(}\sum_{\begin{subarray}{c}a|b\\\
a>A\end{subarray}}\mu(a)\bigg{)}\sideset{}{{}^{\prime}}{\sum}_{\begin{subarray}{c}d\in\mathcal{D}\\\
X/b^{2}\leq|d|\leq(X+h)/b^{2}\end{subarray}}A(b^{2}|d|\sqrt{N},\chi_{b^{2}d}).$

Now we first estimate the partial sum $R$.

######  Proposition 2.1.

$R\ll_{f,\epsilon}A^{-1-\epsilon}(X+h)^{\epsilon+\frac{1}{2}}h^{\frac{1}{2}}+(X+h)^{\frac{1}{2}}hA^{-3+\epsilon}$

###### Proof.

From (1.5), (2.4) and by introducing the factors
$\sum_{d_{1}|b}\mu(d_{1}),\sum_{d_{2}|b}\mu(d_{2})$, as in [4] we have

$\displaystyle A(Q,\chi_{b^{2}d})$
$\displaystyle=\sum_{\begin{subarray}{c}d_{1}|b\end{subarray}}\sum_{\begin{subarray}{c}d_{2}|b\end{subarray}}\mu(d_{1})\mu(d_{2})\frac{\alpha_{d_{1}}\beta_{d_{2}}}{(d_{1}d_{2})^{k}}\chi_{d}(d_{1}d_{2})A(\frac{Q}{d_{1}d_{2}},\chi_{d})$
$\displaystyle=\sum_{\begin{subarray}{c}d_{1}|b\end{subarray}}\sum_{\begin{subarray}{c}d_{2}|b\end{subarray}}\mu(d_{1})\mu(d_{2})\frac{\alpha_{d_{1}}\beta_{d_{2}}}{(d_{1}d_{2})^{k}}\chi_{d}(d_{1}d_{2})\bigg{(}L(k,f,\chi_{d})+O((d_{1}d_{2})^{\frac{1}{2}}|d|Q^{\frac{-1}{2}})\bigg{)}.$

We split the sum into $A,B$, where

$\displaystyle
A=\sum_{\begin{subarray}{c}d_{1}|b\end{subarray}}\sum_{\begin{subarray}{c}d_{2}|b\end{subarray}}\mu(d_{1})\mu(d_{2})\frac{\alpha_{d_{1}}\beta_{d_{2}}}{(d_{1}d_{2})^{k}}\chi_{d}(d_{1}d_{2})L(k,f,\chi_{d})$

and

$\displaystyle
B=\sum_{\begin{subarray}{c}d_{1}|b\end{subarray}}\sum_{\begin{subarray}{c}d_{2}|b\end{subarray}}\mu(d_{1})\mu(d_{2})\frac{\alpha_{d_{1}}\beta_{d_{2}}}{(d_{1}d_{2})^{k}}\chi_{d}(d_{1}d_{2})O\bigg{(}(d_{1}d_{2})^{\frac{1}{2}}|d|Q^{-\frac{1}{2}}\bigg{)}.$

From Deligne’s bound $\alpha_{n},\beta_{n}\leq d(n)n^{(2k-1)/2}$, we have

$\displaystyle
A\leq|L(k,f,\chi_{d})|\sum_{\begin{subarray}{c}d_{1}|b\end{subarray}}\sum_{\begin{subarray}{c}d_{2}|b\end{subarray}}(d_{1}d_{2})^{-\frac{1}{2}+\epsilon}\ll|L(k,f,\chi_{d})|$

and

$\displaystyle
B=O\bigg{(}|d|Q^{\frac{-1}{2}}\sum_{\begin{subarray}{c}d_{1}|b\end{subarray}}\sum_{\begin{subarray}{c}d_{2}|b\end{subarray}}(d_{1}d_{2})^{\frac{\epsilon}{2}}\bigg{)}=O\bigg{(}b^{\epsilon}|d|Q^{\frac{-1}{2}}\bigg{)}.$

Thus

$A(Q,\chi_{b^{2}d})\ll_{f}|L(k,f,\chi_{d})|+b^{\epsilon}|d|Q^{-\frac{1}{2}}.$

Collecting $A$ and $B$ together, and with $Q=b^{2}|d|\sqrt{N}$, $R$ has a
bound

$\displaystyle R$
$\displaystyle\ll_{f}\sum_{(b,4N)=1}\bigg{(}\sum_{\begin{subarray}{c}a|b\\\
a>A\end{subarray}}1\bigg{)}\sideset{}{{}^{\prime}}{\sum}_{\begin{subarray}{c}d\in\mathcal{D}\\\
X/b^{2}\leq|d|\leq(X+h)/b^{2}\end{subarray}}\bigg{(}|L(k,f,\chi_{d})|+b^{-1+\epsilon}|d|^{\frac{1}{2}}\bigg{)}.$

By the Hölder inequality and Theorem 1.2,

$\displaystyle\sideset{}{{}^{\prime}}{\sum}_{\begin{subarray}{c}d\in\mathcal{D}\\\
X/b^{2}\leq|d|\leq(X+h)/b^{2}\end{subarray}}|L(k,f,\chi_{d})|$
$\displaystyle\ll\bigg{(}\sideset{}{{}^{\prime}}{\sum}_{\begin{subarray}{c}d\in\mathcal{D}\\\
X/b^{2}\leq|d|\leq(X+h)/b^{2}\end{subarray}}|L(k,f,\chi_{d})|^{2}\bigg{)}^{\frac{1}{2}}\bigg{(}\sum_{\begin{subarray}{c}X/b^{2}\leq|d|\leq(X+h)/b^{2}\end{subarray}}1\bigg{)}^{\frac{1}{2}}$
$\displaystyle\ll_{\epsilon}((X+h)b^{-2})^{\frac{1+2\epsilon}{2}}h^{\frac{1}{2}}b^{-1}=b^{-2-2\epsilon}(X+h)^{\epsilon+\frac{1}{2}}h^{\frac{1}{2}}$

and

$\displaystyle
b^{-1+\epsilon}\sideset{}{{}^{\prime}}{\sum}_{\begin{subarray}{c}d\in\mathcal{D}\\\
X/b^{2}\leq|d|\leq(X+h)/b^{2}\end{subarray}}|d|^{\frac{1}{2}}\ll
b^{-1+\epsilon}\frac{(X+h)^{\frac{1}{2}}}{b}\frac{h}{b^{2}}=(X+h)^{\frac{1}{2}}hb^{-4+\epsilon}.$

In sum,

$\displaystyle R$
$\displaystyle\ll_{f}\sum_{(b,4N)=1}\bigg{(}\sum_{\begin{subarray}{c}a|b\\\
a>A\end{subarray}}1\bigg{)}\bigg{(}b^{-2-2\epsilon}(X+h)^{\epsilon+\frac{1}{2}}h^{\frac{1}{2}}+(X+h)^{\frac{1}{2}}hb^{-4+\epsilon}\bigg{)}$
$\displaystyle\ll
A^{-1-\epsilon}(X+h)^{\epsilon+\frac{1}{2}}h^{\frac{1}{2}}+(X+h)^{\frac{1}{2}}hA^{-3+\epsilon}.$

∎

We now evaluate $S$. For $(a,4N)=1$ and $d\in\mathcal{D}$, we have

$A(a^{2}|d|\sqrt{N},\chi_{a^{2}d})=\sum_{(n,a)=1}a_{n}n^{-k}\chi_{d}(n)V(\frac{2\pi
n}{a^{2}|d|\sqrt{N}}).$

Write $n=rj^{2}m$, where $r|(4N)^{\infty},(jm,4N)=1$, and $m$ is square-free.
From

$\chi_{d}(n)=\chi_{d}(m)\text{ for }(d,j)=1,$

we obtain

$\displaystyle S$ $\displaystyle=2\sum_{\begin{subarray}{c}a\leq A\\\
(a,4N)=1\end{subarray}}\mu(a)\sum_{\begin{subarray}{c}n=rj^{2}m\\\
(n,a)=1\end{subarray}}a_{n}n^{-k}\sum_{\begin{subarray}{c}(d,j)=1\\\
d\in\mathcal{D}\\\
X/a^{2}\leq|d|\leq(X+h)/a^{2}\end{subarray}}\chi_{d}(m)V(\frac{2\pi
n}{a^{2}|d|\sqrt{N}})$ $\displaystyle=2\sum_{\begin{subarray}{c}a\leq A\\\
(a,4N)=1\end{subarray}}\mu(a)\sum_{\begin{subarray}{c}n=rj^{2}m\\\
(n,a)=1\end{subarray}}a_{n}n^{-k}\sum_{q|j}\mu(q)\sum_{\begin{subarray}{c}dq\in\mathcal{D}\\\
X/a^{2}\leq|d|q\leq(X+h)/a^{2}\end{subarray}}\chi_{dq}(m)V(\frac{2\pi
n}{a^{2}|d|q\sqrt{N}}).$

In the second inequality, we introduced a Möbius factor to relax the
coprimality condition on $d$.

By the Gauss inversion formula,

$\displaystyle S$ $\displaystyle=2\sum_{\begin{subarray}{c}a\leq A\\\
(a,4N)=1\end{subarray}}\mu(a)\sum_{\begin{subarray}{c}n=rj^{2}m\\\
(n,a)=1\end{subarray}}a_{n}n^{-k}\overline{\epsilon_{m}}{m}^{\frac{-1}{2}}\sum_{q|j}\mu(q)$
$\displaystyle\times\sum_{\begin{subarray}{c}dq\in\mathcal{D}\\\
X/a^{2}\leq|d|q\leq(X+h)/a^{2}\end{subarray}}\sum_{2|b|<{m}}\chi_{Nbq}({m})e\bigg{(}\frac{\overline{4N}bd}{m}\bigg{)}V\bigg{(}\frac{2\pi
n}{a^{2}|d|q\sqrt{N}}\bigg{)},$

where

$\displaystyle\overline{\epsilon_{m}}=\begin{cases}i&\text{ if
}m\equiv-1\text{(mod 4)}\\\ 1&\text{ if }m\equiv 1\text{(mod 4)}.\end{cases}$

Set $\Delta=min(\frac{1}{2},a^{2}q(X+h)^{\epsilon-1})$. We split the sum $S$
into

$S=S_{0}+S_{1}+S_{2},$

where the three partial sums are restricted by the conditions
$b=0,0<|b|<\Delta m,\Delta m\leq|b|<m/2$, respectively. Evaluating $|S_{1}|$
and $|S_{2}|$ only requires minor modifications of Iwaniec’s method.

First of all, note that as in [4], $S_{2}\ll_{f}1$.

######  Proposition 2.2.

$S_{1}\ll_{f}A^{2}(X+h)^{\epsilon-\frac{1}{2}}h$

###### Proof.

We have

$\displaystyle S_{1}=$ $\displaystyle 2\sum_{\begin{subarray}{c}a\leq A\\\
(a,4N)=1\end{subarray}}\mu(a)\sum_{\begin{subarray}{c}rj^{2}\\\
(rj,a)=1\end{subarray}}\sum_{q|j}\mu(q)\sum_{\begin{subarray}{c}dq\in\mathcal{D}\\\
X/a^{2}\leq|d|q\leq(X+h)/a^{2}\end{subarray}}$
$\displaystyle\times\sum_{0<|b|}\sum_{\begin{subarray}{c}(m,4Na)=1\\\ m\text{
square -free}\\\
m>|b|/\Delta\end{subarray}}a_{n}n^{-k}\chi_{Nbq}({m})\overline{\epsilon_{m}}{m}^{\frac{-1}{2}}e\bigg{(}\frac{\overline{4N}bd}{m}\bigg{)}V\bigg{(}\frac{2\pi
n}{a^{2}|d|q\sqrt{N}}\bigg{)}.$

Just as in section 8 of [4],

$\displaystyle\sum_{b}\sum_{m}\ll
a^{3}q^{2}r^{-\frac{3}{2}}j^{-3}(X+h)^{\epsilon-\frac{1}{2}}.$

With this estimate, it is straightforward to show

$\displaystyle S_{1}$ $\displaystyle\ll\sum_{\begin{subarray}{c}a\leq A\\\
(a,4N)=1\end{subarray}}\sum_{\begin{subarray}{c}r,j\end{subarray}}\sum_{q|j}\sum_{X/a^{2}q\leq|d|\leq(X+h)/a^{2}q}a^{3}q^{2}r^{-\frac{3}{2}}j^{-3}(X+h)^{\epsilon-\frac{1}{2}}$
$\displaystyle\ll(X+h)^{\epsilon-\frac{1}{2}}h\sum_{\begin{subarray}{c}a\leq
A\end{subarray}}a\sum_{\begin{subarray}{c}r,j\end{subarray}}r^{\frac{-3}{2}}j^{-3}\sum_{q|j}q$
$\displaystyle\ll A^{2}(X+h)^{\epsilon-\frac{1}{2}}h.$

∎

It remains to evaluate $S_{0}$. Because of the condition $b=0$, $S_{0}$ is
written in the following simplified form

$\displaystyle S_{0}$ $\displaystyle=2\sum_{\begin{subarray}{c}a\leq A\\\
(a,4N)=1\end{subarray}}\mu(a)\sum_{\begin{subarray}{c}n=rj^{2}\\\
(n,a)=1\end{subarray}}a_{n}n^{-k}\sum_{q|j}\mu(q)$
$\displaystyle\times\sum_{\begin{subarray}{c}dq\in\mathcal{D}\\\
X/a^{2}q\leq|d|\leq(X+h)/a^{2}q\end{subarray}}V(\frac{2\pi
n}{a^{2}|d|q\sqrt{N}}).$

We then evaluate the innermost sum. We split the inner sum into residue
classes mod $4N$. Each class contributes

$\displaystyle\frac{1}{4N}\int_{X/a^{2}q}^{(X+h)/a^{2}q}V\bigg{(}\frac{2\pi
n}{(a^{2}tq)\sqrt{N}}\bigg{)}dt+O\bigg{(}\bigg{(}1+\frac{n}{X}\bigg{)}^{-\epsilon}\bigg{)}$

by Euler’s summation formula. Hence

$\displaystyle S_{0}$ $\displaystyle=2\sum_{\begin{subarray}{c}a\leq A\\\
(a,4N)=1\end{subarray}}\mu(a)\sum_{\begin{subarray}{c}n=rj^{2}\\\
(n,a)=1\end{subarray}}a_{n}n^{-k}\sum_{q|j}\mu(q)$
$\displaystyle\times\bigg{[}\frac{\gamma(4N)}{4N}\frac{h}{a^{2}q}\int_{0}^{1}V\bigg{(}\frac{2\pi
n}{(ht+X)\sqrt{N}}\bigg{)}dt+O\bigg{(}\bigg{(}1+\frac{n}{X}\bigg{)}^{-\epsilon}\bigg{)}\bigg{]}.$

Here, $\gamma(4N)$ is the order of $\mathcal{D}$. The second term in the inner
sum contributes $O(AX^{1/2+\epsilon})$, by trivial summation over $r,j$.

As in section 9 of [4], the first term in the innermost sum contributes

(2.5)
$\displaystyle\gamma(4N)h\sum_{\begin{subarray}{c}n=rj^{2}\end{subarray}}\frac{a_{n}\phi(j)}{2Nn^{k}j}\sum_{\begin{subarray}{c}a\leq
A\\\
(a,4Nj)=1\end{subarray}}\frac{\mu(a)}{a^{2}}\int_{0}^{1}V\bigg{(}\frac{2\pi
n}{(ht+X)\sqrt{N}}\bigg{)}dt.$

Now, using the identity

$\displaystyle\sum_{\begin{subarray}{c}a\leq A\\\
(a,4Nj)=1\end{subarray}}\frac{\mu(a)}{a^{2}}=\frac{6}{\pi^{2}}\prod_{p\mid
4Nj}(1-p^{-2})^{-1}+O(A^{-1}),$

we split (2.5) into two. The first term is

$\displaystyle\frac{3\gamma(4N)}{\pi^{2}N}\prod_{p|4N}(1-p^{-2})^{-1}h\sum_{\begin{subarray}{c}n=rj^{2}\end{subarray}}\frac{a_{n}\phi(j)}{n^{k}j}\prod_{p|j}(1-p^{-2})^{-1}\int_{0}^{1}V\bigg{(}\frac{2\pi
n}{(ht+X)\sqrt{N}}\bigg{)}dt$
$\displaystyle=\frac{3\gamma(4N)}{\pi^{2}N}\prod_{p|4N}(1-p^{-2})^{-1}h\int_{0}^{1}\bigg{[}\sum_{\begin{subarray}{c}n=rj^{2}\end{subarray}}\frac{a_{n}}{n^{k}}\prod_{p|j}(1+p^{-1})^{-1}V\bigg{(}\frac{2\pi
n}{(ht+X)\sqrt{N}}\bigg{)}\bigg{]}dt.$

Also, the second term contributes $O(A^{-1}hX^{\epsilon})$, by trivial
summation over $r,j$.

Thus, finally,

$\displaystyle S_{0}$
$\displaystyle=\frac{3\gamma(4N)}{\pi^{2}N}\prod_{p|4N}(1-p^{-2})^{-1}h\int_{0}^{1}\mathcal{B}(ht+X)dt$
$\displaystyle+O(A^{-1}hX^{\epsilon}+AX^{1/2+\epsilon}),$

where

$\displaystyle\mathcal{B}(x)=\sum_{\begin{subarray}{c}n=rj^{2}\end{subarray}}\frac{a_{n}}{n^{k}}\prod_{p|j}(1+p^{-1})^{-1}V\bigg{(}\frac{2\pi
n}{x}\bigg{)}.$

By shifting the line of integration to $(-\frac{1}{5})$,

$\displaystyle\mathcal{B}(x)$ $\displaystyle=\frac{1}{2\pi
i}\int_{(4/5)}\frac{\Gamma(k+s)}{\Gamma(k)}\sum_{\begin{subarray}{c}n=rj^{2}\end{subarray}}\frac{a_{n}}{n^{k+s}}\prod_{p|j}(1+p^{-1})^{-1}\bigg{(}\frac{2\pi}{x}\bigg{)}^{-s}\frac{ds}{s}$
$\displaystyle=-\frac{1}{2\pi
i}\int_{(-1/5)}\frac{\Gamma(k+s)}{\Gamma(k)}\bigg{(}\frac{2\pi}{x}\bigg{)}^{-s}L_{f}(k+s)\frac{ds}{s}+L_{f}(k)$

Note that $L_{f}(s)$ appeared in [4] as $L(s)$ for $k=1$, and it was later
generalized in [8] as $L_{f,1}(s)$ for all positive integer $k$. Recall that
$L_{f}(k)\neq 0$ and $L_{f}(s)$ is polynomially bounded if $\Re(s)>k-1/4$.
Thus we have

(2.6) $\displaystyle\mathcal{B}(x)=L_{f}(k)+O_{f}(x^{-\frac{1}{5}}).$

From (2.6), we obtain

$\displaystyle
S_{0}=\frac{3\gamma(4N)}{\pi^{2}N}\prod_{p|4N}(1-p^{-2})^{-1}L_{f}(k)h+O(h^{\frac{4}{5}}+A^{-1}hX^{\epsilon}+AX^{\frac{1}{2}+\epsilon}).$

Now we collect all the error terms of $S_{f}(X,h)$. They are

$\displaystyle O_{f,\epsilon}(h^{\frac{4}{5}}$
$\displaystyle+A^{-1}hX^{\epsilon}+AX^{\frac{1}{2}+\epsilon}+A^{2}(X+h)^{\epsilon-\frac{1}{2}}h$
$\displaystyle+A^{-1-\epsilon}(X+h)^{\epsilon+\frac{1}{2}}h^{\frac{1}{2}}+(X+h)^{\frac{1}{2}}hA^{-3+\epsilon}).$

Let $A=X^{\frac{1}{4}-\epsilon}$. Then for $X^{3/4+\epsilon}\leq h\leq X$, the
main term dominates the error terms. We have

$\displaystyle S_{f}(X,h)$
$\displaystyle=\frac{3\gamma(4N)}{\pi^{2}N}\prod_{p|4N}(1-p^{-2})^{-1}L_{f}(k)h+O_{f,\epsilon}(hX^{-\epsilon}).$

## Acknowledgement

This work was supported by the National Research Foundation of Korea (NRF)
grant funded by the Korea government(MSIT)(No. 2019R1A2C108860913, No.
2020R1A4A1016649). The author thanks his advisor Hae-Sang Sun for helpful
comments and discussions. The author also thanks to Yoonbok Lee, Jaeseong Kwon
and Seongjae Han for their helpful comments and discussions.

## References

  * [1] Antal Balog and Ken Ono.  The chebotarev density theorem in short intervals and some questions of serre.  Journal of Number Theory, 91(2):356–371, 2001.
  * [2] Soumya Das and Satadal Ganguly.  Gaps between nonzero fourier coefficients of cusp forms.  Proceedings of the American Mathematical Society, 142(11):3747–3755, 2014.
  * [3] Kazuyuki Hatada.  Eigenvalues of hecke operators on $SL(2,\mathbb{Z})$.  Mathematische Annalen, 239(1):75–96, 1979.
  * [4] Henryk Iwaniec.  On the order of vanishing of modular $L$-functions at the critical point.  Journal de théorie des nombres de Bordeaux, 2(2):365–376, 1990\.
  * [5] Benjamin Justus.  On the moments of central values of modular $L$-functions.  PhD thesis, Concordia University, 2008.
  * [6] Emmanuel Kowalski, Olivier Robert, and Jie Wu.  Small gaps in coefficients of $L$-functions and $\mathfrak{B}$-free numbers in short intervals.  Revista Matemática Iberoamericana, 23(1):281–326, 2007.
  * [7] Narasimha Kumar.  On the gaps between non-zero fourier coefficients of cusp forms of higher weight.  The Ramanujan Journal, 45(1):95–109, 2018.
  * [8] Wenzhi Luo and Dinakar Ramakrishnan.  Determination of modular forms by twists of critical $L$-values.  Inventiones mathematicae, 130(2):371–398, 1997.
  * [9] M Ram Murty and V Kumar Murty.  Mean values of derivatives of modular $L$-series.  Annals of Mathematics, 133(3):447–475, 1991.
  * [10] Ken Ono and Christopher Skinner.  Non-vanishing of quadratic twists of modular $L$-functions.  Inventiones mathematicae, 134(3):651–660, 1998.
  * [11] Ariel Pacetti.  Shimura correspondence and central values of twisted $L$-series.  Proceedings of the “Segundas Jornadas de Teorıa de Números”(Madrid, 2007), 1:21, 2008.
  * [12] Alberto Perelli and J Pomykała.  Averages of twisted elliptic $L$-functions.  Acta Arithmetica, 80:149–163, 1997.
  * [13] Jacek Pomykała.  Non-vanishing of $n$-th derivatives of twisted elliptic $L$-functions in the critical point.  Journal de théorie des nombres de Bordeaux, 9(1):1–10, 1997\.
  * [14] Jean-Pierre Serre.  Quelques applications du théoreme de densité de chebotarev.  Publications Mathématiques de l’IHÉS, 54:123–201, 1981\.


# Robust Expected Information Gain for Optimal Bayesian Experimental Design
Using Ambiguity Sets

Jinwoo Go Computational Science and Engineering Dept.
Georgia Institute of Technology
Atlanta, Georgia, USA  Tobin Isaac  Computational Science and Engineering
Dept.
Georgia Institute of Technology
Atlanta, Georgia, USA

###### Abstract

The ranking of experiments by expected information gain (EIG) in Bayesian
experimental design is sensitive to changes in the model’s prior distribution,
and the approximation of EIG yielded by sampling will have errors similar to
the use of a perturbed prior. We define and analyze _robust expected
information gain_ (REIG), a modification of the objective in EIG maximization
by minimizing an affine relaxation of EIG over an ambiguity set of
distributions that are close to the original prior in KL-divergence. We show
that, when combined with a sampling-based approach to estimating EIG, REIG
corresponds to a ‘log-sum-exp’ stabilization of the samples used to estimate
EIG, meaning that it can be efficiently implemented in practice. Numerical
tests combining REIG with variational nested Monte Carlo (VNMC), adaptive
contrastive estimation (ACE) and mutual information neural estimation (MINE)
suggest that in practice REIG also compensates for the variability of under-
sampled estimators.

##  1 Introduction

Bayesian Experimental Design (BED) is a probabilistic framework for selecting
experiments to learn about one or more uncertain variables. Within BED, the
most popular criterion for ranking experiments is by Expected Information Gain
(EIG), which estimates from current knowledge, encoded in a prior
distribution, how informative a particular experiment is likely to be. This
framework is used in diverse applications across many disciplines [Ryan et
al., 2016a], in natural sciences [Huan, 2010], social sciences [Embretson and
Reise, 2013], and in machine learning and data analysis [Foster et al., 2020].

The sensitivity of experimental design to misspecification of the prior
distribution has been described in [DasGupta and Studden, 1991, Ryan et al.,
2016b], even in settings where the information gain is computable in closed
form. For more complex models, EIG can only be estimated numerically, which
may further affect the reliability of the computed rankings of experiments
(or, in the case of continuously parameterized experiments, the gradient of
EIG). EIG is by definition an expectation of an expectation, so general
purpose estimates, such as Nested Monte Carlo (NMC) estimation [Ryan, 2003],
can be expensive, slow to converge, and sensitive to underconverged sample
estimates.

To address the issues above that affect the reliability of EIG estimates in
BED, we introduce a quantity we call _robust expected information gain_ (REIG)
as a probability-theoretic way of ranking experiments by their expected
information gain for some worst-case small perturbation of the prior. We also
show through convex analysis that the estimation of REIG is a simple post-
processing of the samples generated by an NMC-like estimator. As a result, our
methodology is applicable with many existing methods, which we demonstrate in
section 7 by applying REIG to samples generated by three recent popular EIG
estimators [Foster et al., 2020, 2019, Kleinegesse and Gutmann, 2020].

##  2 Background and Notation

We use $\theta\in\Theta$ to indicate a choice of parameters for a model from a
set of possible parameters, and we assume a reference prior probability
distribution of $\theta$ which has a measurable density function $p(\theta)$,
so that we can write $E_{p(\theta)}[f]=\int_{\Theta}f(\theta)p(\theta)\
d\theta.$ We let $\xi\in\Xi$ be a potential experiment from a class of
experiments, which has an outcome variable $y(\xi)$. The experiment $\xi$ is
modeled by the likelihood function $p(y|\theta,\xi)$, which for each choice of
$\theta$ defines a measurable probability density function of $y$. Our
interest is in models where the densities $p(\theta)$ and $p(y|\theta,\xi)$
can be efficiently computed, and where samples can be drawn from $p(\theta)$
and from $p(y|\theta,\xi)$ for each $(\theta,\xi)$, so that the joint prior
distribution $p(\theta,y|\xi)=p(\theta)p(y|\theta,\xi)$ also has a computable
density function and can be sampled. We use the notation
$p(y|\xi)=E_{p(\theta)}[p(y|\theta,\xi)]$ for the marginal distribution of the
outcome $y$.

###  2.1 Prior Uncertainty

While it seems recursive to consider uncertainty in the prior distribution
used in Bayesian inference, the prior distribution is in many settings not
determined from first principles or an existing population of data. In these
cases the choice of prior is often dictated by what is required to make a
computation tractable or simple, by invariance principles, or by an attempt to
be noninformative [Stark, 2015]. But notions of noninformative priors do not
scale to high dimensions [Yang and Berger, ], and even in low dimensions
priors that are close under a weak topology like total variation can have
diverging posteriors for the same observations [Owhadi et al., 2015].

So in this work we will consider sets of prior distributions $q(\theta)$ other
than the reference $p(\theta)$, but we only consider $q(\theta)$ that are
absolutely continuous with respect to $p(\theta)$. Although methods similar to
ours are used to handle model uncertainty [Shapiro et al., 2021], we treat
$p(y|\theta,\xi)$ as certain: the only uncertainty we consider is in the prior
distribution of $\theta$. When we extend the notation of derived distributions
from $p(\theta)$ to another $q(\theta)$, then, it is always with the same
likelihood: the joint prior $q(\theta,y|\xi)=q(\theta)p(y|\theta,\xi)$, the
marginal $q(y|\xi)=E_{q(\theta)}[p(y|\theta,\xi)]$, etc.

###  2.2 Expected Information Gain

We use $D_{\mathrm{KL}}(p(\theta)\|q(\theta))$ to denote the Kullback-Leibler
divergence,
$D_{\mathrm{KL}}(p(\theta)\|q(\theta))=\int_{\Theta}p(\theta)\log(p(\theta)/q(\theta))\
d\theta$. The expected information gain of experiment $\xi$ is defined to be
the expectation over the marginal distribution of outcomes $p(y|\xi)$ of the
KL-divergence from the Bayesian posterior distribution $p(\theta|y,\xi)$ to
the prior $p(\theta)$. Because the prior is not fixed in this work, we
consider EIG to be a function of both the prior $p(\theta)$ and the experiment
$\xi$,

$I(p,\xi)=E_{p(y|\xi)}[D_{\mathrm{KL}}(p(\theta|y,\xi)\|p(\theta))].$ (1)

Bayesian optimal experimental design seeks the experiment $\xi^{*}$ that
maximizes this quantity,

$\xi^{*}=\arg\max_{\xi\in\Xi}\ I(p,\xi).$

Although the form of EIG in (1) is the most intuitive for Bayesian
experimental design, other equivalent definitions map more directly on the
robust variant we introduce in section 5 and the sampling-based estimators in
section 6. The EIG of an experiment is also the mutual information between
$\theta$ and $y$,

$\displaystyle I(p,\xi)$
$\displaystyle=D_{\mathrm{KL}}(p(\theta,y|\xi)\|p(\theta)p(y|\xi))$ (2)
$\displaystyle=E_{p(\theta)}[D_{\mathrm{KL}}(p(y|\theta,\xi)\|p(y|\xi))]$ (3)
$\displaystyle=E_{p(\theta,y|\xi)}\big{[}\log\frac{p(y|\theta,\xi)}{p(y|\xi)}\big{]}.$
(4)

The last form in (4) is the preferred form for many methods that estimate EIG
when the likelihood $p(y|\theta,\xi)$ and prior $p(\theta)$ can be evaluated
directly. We discuss methods for estimating $I(p,\xi)$ and other related
quantities in section 6.

##  3 A Simple Example with Two Experiments

Suppose a doctor has two blood tests for Condition X: test A has a
$10^{-14}\%$ chance of a false negative but a $50\%$ chance of a false
positive, and test B has an $\approx 18.4\%$ chance of a false negative and
the same chance of a false positive. If the doctor estimates the prior
probability that a patient has condition X is $50\%$, it turns out that both
tests have the same expected information gain of $\approx 0.22$ nats. If that
prior probability could be mistaken, however, the two tests have different
EIGs for prior probabilities in the vicinity of $50\%$, shown in fig. 1. If
the prior probability the patient has Condition X is actually $>50\%$, then
test A has a greater EIG than test B, and vice versa if it is $<50\%$.

Figure 1: The expected information gain of two tests for condition X depends
on the condition’s prior probability.

We interpret these results as follows: For test A, if the prior probability of
Condition X is $>50\%$, then a negative test result is surprising because
there are essentially no false negatives, but if the prior probability is
$<50\%$, then a positive result is less surprising because it has a high
probability of being a false positive. So in comparison to test B, the EIG of
test A is more sensitive to the choice of prior. In any neighborhood of
$p=50\%$, there are priors where test A is expected to be less informative
than test B. So one could argue that a risk-averse doctor, who would maximize
how informative the test would be in the worst case, should select test B.

##  4 Ambiguity Sets

In the simple example above, we used a range of prior probabilities for the
model parameter to argue that some tests are less locally sensitive to
perturbations of the prior distribution. To generalize this idea from a simple
discrete example to other probability distributions, we rely on the notion of
an _ambiguity set_ [Bayraksan and Love, 2015], which is a set of distributions
that are not far from a reference prior $p(\theta)$ in some statistical
distance. We use KL-divergence as distance, so our ambiguity set with radius
$\epsilon$ centered at reference distribution $p(\theta)$ is

$\mathcal{A}(\epsilon,p)=\\{q:D_{\mathrm{KL}}(q(\theta)\|p(\theta))\leq\epsilon\\}.$

KL-divergence as a distance works well with Bayesian optimal experimental
design, because the KL-divergence appears in the definition of EIG, and
because the set $\mathcal{A}(\epsilon,p)$ is defined as a convex subset of
positive measurable functions $q(\theta)$ with just two conditions:
$\int_{\Theta}q(\theta)\ d\theta=1$ and
$D_{\mathrm{KL}}(q(\theta)\|p(\theta))\leq\epsilon$. Thus the minimization
over $q\in\mathcal{A}(\epsilon,p)$ of a well-behaved convex function $f(q)$,
which appears difficult because $\mathcal{A}(\epsilon,p)$ of the infinite
dimensionality of space of measurable functions, transforms into an equivalent
dual convex program with only two variables.

We direct the interested reader to [Shapiro, 2017] for additional details:
here we summarize the results that are important for this work. If the
objective function of interest $f(q)$ is the expectation under $q(\theta)$ of
a measurable quantity of interest $Z(\theta)$, then it is an affine function
of $q$ and we may use duality to simplify the maximization or minimization of
$E_{q(\theta)}[Z(\theta)]$ over $\mathcal{A}(\epsilon,p(\theta))$ into a dual
problem with only one variable $\lambda\geq 0$. The maximization problem

$R_{\epsilon}=\sup_{q\in\mathcal{A}(\epsilon,p)}[Z(\theta)]$ (5)

can be solved in dual form as

$R_{\epsilon}=\inf_{\lambda\geq 0}\lambda\epsilon+\lambda\log
E_{p(\theta)}[\exp(\lambda^{-1}Z(\theta))],$ (6)

and the minimization problem

$M_{\epsilon}=\inf_{q\in\mathcal{A}(\epsilon,p)}[Z(\theta)]$ (7)

can be solved in dual form as

$M_{\epsilon}=-\inf_{\lambda\geq 0}\lambda\epsilon+\lambda\log
E_{p(\theta)}[\exp(-\lambda^{-1}Z(\theta))].$ (8)

It is important to note that these are non-parametric results. Given a
parameterized family of priors $\\{p(\theta;\psi)\\}_{\psi}$, the gradient
$\nabla_{\psi}E_{p(\theta;\psi)}[Z(\theta)]$ is sufficient to compute the
optimizer in $\mathcal{A}(\epsilon,p)$ of (5) or (7) within the parametric
family because the objective is affine. But (6) and (8) allow us to compute
the optimal objective value over the entire ambiguity set without explicitly
computing an optimal distribution in the set.

##  5 Robust Bayesian Experimental Design

The insight of the example of section 3 was that a risk-averse approach to
experimental design that allows for some uncertainty in the prior distribution
would select the experiment that maximizes the worst-case EIG in the vicinity
of the reference prior $p(\theta)$. Using the ambiguity set
$\mathcal{A}(\epsilon,p(\theta))$ from section 4 to define the vicinity of
$p(\theta)$, we first formalize the worst-case EIG as
$I_{\epsilon}^{\mathrm{true}}(p,\xi)$, the _true robust expected information
gain with radius $\epsilon$_,

$I_{\epsilon}^{\mathrm{true}}(p,\xi)=\inf_{q\in\mathcal{A}(\epsilon,p)}I(q,\xi).$
(9)

The experiment that maximizes this quantity is

$\displaystyle\xi^{*}_{\mathrm{REIG},\epsilon,\mathrm{true}}$
$\displaystyle=\arg\max_{\xi\in\Xi}\ I_{\epsilon}^{\mathrm{true}}(p,\xi)$
$\displaystyle=\arg\max_{\xi\in\Xi}\inf_{q\in\mathcal{A}(\epsilon,p)}I(q,\xi).$

This optimization problem in this definition has a clear meaning, but we note
that $I(q,\xi)$ is a quantity that is concave in $q$, so it can have multiple
local minima in the convex ambiguity set $\mathcal{A}(\epsilon,p)$ and the
duality framework of section 4 cannot be applied directly.

###  5.1 Affine Expected Information Gain Approximation

To define a relaxation to a tractable problem, we split $I(q,\xi)$ into two
contributions, one with the marginal distribution of $y$ fixed by the
reference prior, $y\sim p(y|\xi)$, and the other a correction that is the
divergence between $p(y|\xi)$ and $q(y|\xi)$,

$\displaystyle I(q,\xi)$
$\displaystyle=E_{q(\theta,y|\xi)}\big{[}\log\big{(}\frac{p(y|\theta,\xi)}{p(y|\xi)}\frac{p(y|\xi)}{q(y|\xi)}\big{)}\big{]}$
(10)
$\displaystyle=E_{q(\theta,y|\xi)}\big{[}\log\frac{p(y|\theta,\xi)}{p(y|\xi)}\big{]}$
$\displaystyle\phantom{=}{}-D_{\mathrm{KL}}(q(y|\xi)\|p(y|\xi)).$

We denote the first term in this difference

$I_{\mathrm{aff}}(q,\xi;p)=E_{q(\theta,y|\xi)}\big{[}\log\frac{p(y|\theta,\xi)}{p(y|\xi)}\big{]},$
(11)

because it is an approximation to $I(q,\xi)$ that is affine and exact when
$q=p$. By the concavity of EIG with respect to its first argument q,

$I_{\mathrm{aff}}(q,\xi;p)\geq I(q,\xi)\quad\text{for all }q.$ (12)

Due to the data processing inequality that the mutual information between two
random variables cannot increase by a deterministic or random transformation
of the arguments, the error in $I(q,\xi;p)$ is bounded by

$\displaystyle|I(q,\xi)-I_{\mathrm{aff}}(q,\xi;p)|$ $\displaystyle\leq
D_{\mathrm{KL}}(q(y|\xi)\|p(y|\xi))$ (13) $\displaystyle\leq
D_{\mathrm{KL}}(q(\theta)\|p(\theta)).$

In fact $I_{\mathrm{aff}}(q,\xi;p)$ is the affine approximation to $I(q,p)$
that is tangent at $q=p$.

######  Theorem 1.

The function $I_{\mathrm{aff}}(q,\xi;p)$ from (11) is tangent to $I(q,\xi)$ at
$q=p$ for every design $\xi$.

###### Proof.

It is sufficient to show that the difference between the two functions, which
by (10) is $D_{\mathrm{KL}}(q(y|\xi)\|p(y|\xi))$, is gradient free at $q=p$.

We first calculate the gradient with respect to $\xi$: by the chain rule
applied to (4), the derivative in the direction $\hat{\xi}$ is

$\displaystyle\nabla_{\xi}D_{\mathrm{KL}}(q(y|\xi)\|p(y|\xi))[\hat{\xi}]$
$\displaystyle={}$
$\displaystyle\nabla_{q(y|\xi)}D_{\mathrm{KL}}(q(y|\xi)\|p(y|\xi))[\nabla_{\xi}q(y|\xi)[\hat{\xi}]]$
$\displaystyle{}-E_{q(y|\xi)}\big{[}\frac{\nabla_{\xi}p(y|\xi)[\hat{\xi}]}{p(y|\xi)}\big{]}.$

For general distributions $Q$ and $P$ the KL divergence satisfies
$\nabla|_{Q=P}D_{\mathrm{KL}}(Q\|P)=0$, so the first term vanishes when $q=p$.
In the second term, when $q=p$ the denominator cancels with the measure and we
have

$\displaystyle
E_{q(y|\xi)}\big{[}\frac{\nabla_{\xi}p(y|\xi)[\hat{\xi}]}{p(y|\xi)}\big{]}|_{q=p}$
$\displaystyle={\textstyle\int}\nabla_{\xi}p(y|\xi)[\hat{\xi}]\ dy$
$\displaystyle=\nabla_{\xi}(E_{p(y|\xi)}[1])[\hat{\xi}]=0,$

where we use the fact that $E_{p(y|\xi)}[1]=1$ for all $\xi$.

Finally, we can see that the gradient with respect to $q(\theta)$ in the
direction $\hat{q}(\theta)$ is

$\displaystyle\nabla_{q(\theta)}D_{\mathrm{KL}}(q(y|\xi)\|p(y|\xi))[\hat{q}(\theta)]$
$\displaystyle={}$
$\displaystyle\nabla_{q(y|\xi)}D_{\mathrm{KL}}(q(y|\xi)\|p(y|\xi))[\nabla_{q(\theta)}q(y|\xi)[\hat{q}(\theta)]],$

which vanishes at $q=p$ for the same reasons as above. ∎

###  5.2 Robust Expected Information Gain (REIG)

Having shown that $I_{\mathrm{aff}}(q,\xi;p)$ is a good approximation to
$I(q,\xi)$ near the reference prior $p(\theta)$, we now use it to define a
robust quantity that approximates $I_{\epsilon}^{\mathrm{true}}$, which we
refer to simply as $I_{\epsilon}$,

$I_{\epsilon}(p,\xi)=\inf_{q\in\mathcal{A}(\epsilon,p)}I_{\mathrm{aff}}(q,\xi;p).$
(14)

By the properties established in (12), (13), and theorem 1, we have the
following relationships between $I_{\epsilon}^{\mathrm{true}}$ and
$I_{\epsilon}$:

$\displaystyle I(p,\xi)\geq I_{\epsilon}(p,\xi)\geq
I_{\epsilon}^{\mathrm{true}}(p,\xi);$ (15)
$\displaystyle|I_{\epsilon}(p,\xi)-I_{\epsilon}^{\mathrm{true}}(p,\xi)|\leq\epsilon;$
(16)
$\displaystyle|I_{\epsilon}(p,\xi)-I_{\epsilon}^{\mathrm{true}}(p,\xi)|\in
O(\epsilon^{2}).$ (17)

These facts suggest an experiment that maximizes $I_{\epsilon}(p,\xi)$,

$\xi^{*}_{\mathrm{REIG},\epsilon}=\arg\max_{\xi\in\Xi}I_{\epsilon}(p,\xi),$

has similar robustness to $\xi^{*}_{\mathrm{REIG},\epsilon,\mathrm{true}}$
over perturbations of the the prior $p(\theta)$, as long as the radius
$\epsilon$ of the ambiguity set is not too large.

###  5.3 Computation of $I_{\epsilon}(p,\xi)$ via duality

We have selected $I_{\epsilon}(p,\xi)$ as our robust quantity to optimize
because (14) can be optimized by the dual transformation described in section
4. Applying (8) to EIG in the form (3), we have $I_{\epsilon}(p,\xi)=$

$-\inf_{\lambda\geq 0}\lambda\epsilon+\lambda\log
E_{p(\theta)}\big{[}\exp\big{(}-\frac{D_{\mathrm{KL}}(p(y|\theta,\xi)\|p(y|\xi))}{\lambda}\big{)}\big{]}.$
(18)

We will show in section 6 that this 1D convex optimization problem can be
solved efficiently by a small adaptation of existing EIG estimators.

###  5.4 Related Design Criteria

Our definition of $I_{\epsilon}$ was motivated by a risk-aversion argument in
favor of the design with the best worst-case EIG in a neighborhood. Because
the approximation $I_{\mathrm{aff}}(q,\xi;d)$ is affine, however, maximization
over the ambiguity set can also be solved by duality. This means that the same
methodology can be used to define a risk-loving strategy for experimental
design, which selects the experiment that has the highest EIG for some prior
in the ambiguity set. We call this criterion $I_{\epsilon,\max}(p,\xi)=$

$\inf_{\lambda\geq 0}\lambda\epsilon+\lambda\log
E_{p(\theta)}\big{[}\exp\big{(}\frac{D_{\mathrm{KL}}(p(y|\theta,\xi)\|p(y|\xi))}{\lambda}\big{)}\big{]}.$
(19)

This criterion is used in section 7 to counteract biased underestimation of
EIG by some estimators.

Last, we note that our decision to limit the uncertainty in the models of the
experiments to just the prior $p(\theta)$ and not the likelihood
$p(y|\theta,\xi)$ is arbitrary, at least from the perspective of the methods
we have developed. An ambiguity set $\mathcal{A}(\epsilon,p(\theta,y|\xi))$
can be centered around the joint prior of the model $p(\theta,y|\xi)$, and an
affine approximation can be taken that would allow for optimization over that
ambiguity set via duality. The result would be an even more conservative
quantity, $I_{\epsilon,\text{joint}}(p,\xi)=$

$-\inf_{\lambda\geq 0}\lambda\epsilon+\lambda\log
E_{p(\theta,y|\xi)}\big{[}\exp\big{(}\lambda^{-1}\log\frac{p(y|\xi)}{\log
p(y|\theta,\xi)}\big{)}\big{]}.$ (20)

We will not explore this criterion more in this work.

##  6 REIG Estimation via Sampling

The design of efficient estimators for EIG has been the subject of much
research, in part because their use in experimental design is computationally
demanding. The combination of nested iterations to estimate the densities of
implicitly defined distributions, to evaluate the expectation of the EIG, and
finally to optimize that quantity lead to many passes over the problem data as
well as many model evaluations.

When introducing an implicitly defined quantity like $I_{\epsilon}$, we should
be leery of adding another nested loop to the calculation. This is why we
immediately discounted the design criterion $I_{\epsilon}^{\mathrm{true}}$
from (9), which would require optimization in the original variables
parameterizing $p(\theta)$, which could be numerous.

###  6.1 Constructing a REIG Estimator

When both the prior $p(\theta)$ and the likelihood $p(y|\theta,\xi)$ can be
sampled directly, sampling-based approaches to estimating EIG often have a
two-level structure: an inner estimator is defined for a fixed $\theta$ and/or
$y$ in the integrated quantity — either
$D_{\mathrm{KL}}(p(\theta|y,\xi)\|p(\theta))$ in (1),
$D_{\mathrm{KL}}(p(y|\theta,\xi)\|p(y|\xi))$ in (3), or
$\log(p(y|\theta,\xi)/p(y|\xi)$ in (4) — and an outer Monte Carlo estimator
over either $p(\theta)$ or $p(\theta,y|\xi)$ calls the inner estimator for
each generated $\theta$ or $(\theta,y)$.

This basic paradigm maps closely onto the dual formulation of $I_{\epsilon}$
in (18), in a method we sketch in algorithm 1 that defines a REIG estimator
$\hat{I}_{\epsilon}$.

  1. 1.

Draw $N_{1}$ i.i.d. samples $\\{\theta_{i}\\}_{i=1}^{N_{1}}$ from $p(\theta)$.

  2. 2.

For each $\theta_{i}$, use estimator $\tilde{D}(\theta,\xi)$ to compute an
estimate $d_{i}\leftarrow\tilde{D}(\theta_{i},\xi)$ of
$D_{\mathrm{KL}}(p(y|\theta,\xi)\|p(y|\xi))$.

  3. 3.

Solve the 1D convex optimization problem

$M_{\epsilon}=\inf_{\lambda\geq
0}\lambda\epsilon+\lambda\log\frac{1}{N_{1}}\sum_{i=1}^{N_{1}}\exp(-\lambda^{-1}d_{i})$
(21)

and return $-M_{\epsilon}$.

Algorithm 1 $I_{\epsilon}$ Estimation via Sampling

In this approach the inner estimator is called $N_{1}$ times in step 2, which
is the same number of times it would have been called to compute the EIG
estimator $\frac{1}{N_{1}}\sum_{i=1}^{N_{1}}d_{i}$, but those estimates are
saved and treated as an empirical distribution, so that the optimization
problem in step 3 solves (18) by sample average approximation (SAA) instead of
stochastic approximation (SA). The assumption is that this 1D convex problem
will be solved quickly and the dominant cost in algorithm 1 is the cost of
computing the KL-divergence estimators
$d_{i}\leftarrow\tilde{D}(\theta_{i},\xi)$.

Although the inner optimization in step 3 is solved by SAA, we note that
algorithm 1 can be used within either SAA or SA for the optimization over
$\xi$, depending on whether the samples in step 1 are reused or not. The
derivative $\nabla_{\xi}\hat{I}_{\epsilon}(p,\xi)$ can be computed as
$-\nabla_{d}M_{\epsilon}\cdot\nabla_{\xi}d$, where $d$ is the vector of
$d_{i}$ estimates from step 2. The partial derivatives $\nabla_{\xi}d$ are
also present in computing the gradients of EIG estimators, so existing methods
for this term can be reused: $\nabla_{d}M_{\epsilon}$ are the only additional
derivatives needed for REIG. Letting $\lambda^{*}$ be the optimal value and
letting $L(d)=\log\frac{1}{N_{1}}\sum_{i=1}^{N_{1}}\exp(-d_{i})$, if
$\lambda^{*}\neq 0$ then $\nabla_{d}M_{\epsilon}=\nabla_{d}L(d/\lambda^{*})$.
If $\lambda^{i}=0$ there there is some $i^{*}=\arg\min_{i}d_{i}$ and
$M_{\epsilon}=-d_{i}$, so that either $-e_{i}=\nabla_{d}M_{\epsilon}$ or
$-e_{i}\in\partial_{d}M_{\epsilon}$ if $i^{*}$ is not unique.

###  6.2 EIG Estimators

There are many possible choices for the estimator $\tilde{D}(\theta,\xi)$ of
$D_{\mathrm{KL}}(p(y|\theta,\xi)\|p(y|\xi))$. In fact, any EIG estimator
$\hat{D}(p(\theta),p(y|\theta,\xi))$ that accepts general prior distributions
can be used by running a separate instance of $\hat{D}$ for each $\theta_{i}$
with the prior distribution $\theta\sim\delta_{\theta_{i}}$. In practice, this
approach would have poor performance because sequestering the samples
$\theta_{i}$ into separate estimator instances would not allow for
vectorization across samples. Vectorization and batching are best exploited in
a nested EIG estimator if all instances of the inner estimator have the same
hyperparameters to maximize throughput.

In the experiments in section 7, we have primarily used the estimators
developed in [Foster et al., 2019, 2020, Kleinegesse and Gutmann, 2020]. The
first is the _Variational Nested Monte Carlo_ (VNMC) estimator, which is based
on the form of EIG in (4). It draws $N$ samples $(\theta_{i},y_{i})\sim
p(\theta,y|\xi)$ from the joint prior, evaluates $\log
p(y_{i}|\theta_{i},\xi)$ directly, and then uses an inner estimator for $\log
p(y_{i}|\xi)$. That estimator is based on the identity

$\log p(y|\xi)=\log
E_{q(\theta)}\big{[}\frac{p(\theta)}{q(\theta|y,\xi)}p(y|\theta,\xi)\big{]},$

where $q(\theta|y,\xi)$ can be any distribution that is absolutely continuous
with respect to $p(\theta)$, but the variance of that expectation is lower the
closer $q(\theta|y,\xi)$ is to the posterior distribution $p(\theta|y,\xi)$.
$M$ samples $\\{\theta_{i}^{j}\\}_{j=1}^{M}$ are drawn from $q(\theta|y,\xi)$
for each $y_{i}$, resulting in the EIG estimator

$\hat{I}^{\text{VNMC}}=\sum_{i=1}^{N}\log\frac{p(y_{i}|\theta_{i},\xi)}{\frac{1}{M}\sum_{j=1}^{M}\frac{p(\theta_{i}^{j})}{q(\theta_{i}^{j}|y_{i},\xi)}p(y_{i}|\theta_{i}^{j},\xi).}$
(22)

This estimator is consistent in the limit as $M\to\infty$, but for finite $M$
is in expectation an upper bound for EIG. For details see Foster et al.
[2019].

The second estimator that we use is the _Adaptive Contrastive Estimator_ (ACE)
from [Foster et al., 2020]. In description it is almost identical to VNMC,
except that in estimating $p(y_{i}|\xi)$ we add to the samples
$\\{\theta_{i}^{j}\\}_{j=1}^{M}$ the original sample
$\theta_{i}^{0}=\theta_{i}$ drawn from $p(\theta)$ that generated $y_{i}$. The
result is

$\hat{I}^{\text{ACE}}=\sum_{i=1}^{N}\log\frac{p(y_{i}|\theta_{i},\xi)}{\frac{1}{M+1}\sum_{j=0}^{M}\frac{p(\theta_{i}^{j})}{q(\theta_{i}^{j}|y_{i},\xi)}p(y_{i}|\theta_{i}^{j},\xi).}$
(23)

This is also a consistent estimator of EIG, but the addition of the prior-
samples $\theta_{i}^{0}$ to the estimator for $p(y_{i}|\xi)$ makes it in
expectation a lower bound for EIG for finite $M$: see [Foster et al., 2020]
for more details.

The last estimator that we use is the _Mutual Information Neural Estimation_
(MINE), which trains the ratio, $\frac{p(y|\theta,\xi}{p(y|\xi)}$, with
samples and estimate the EIG with SAA method [Kleinegesse and Gutmann, 2020],

$\hat{I}^{\text{MINE}}=\sum_{i}^{N}[T_{\psi}(\theta_{i},y_{i})-e^{T_{\psi}(\theta_{i},y_{i}^{*})-1}],$
(24)

where $y_{i}^{*}$ represents the shuffled $y$ and $T$ is a neural network with
parameters $\psi$.

The outer loop for ACE, VNMC and MINE methods draw from $p(\theta,y|\xi)$ or
$p(\theta)p(y|\xi)$, and in the implementations provided by the authors one
sample from $p(y|\theta_{i},\xi)$ is drawn for each of $N$ samples
$\theta_{i}$ drawn from $p(\theta)$. To adapt these methods to the needs of
our algorithm 1, we split $N$ into $N=N_{1}N_{2}$, and draw $N_{2}$ samples
from $p(y|\theta_{i},\xi)$ for each of $N_{1}$ samples drawn from $p(\theta)$.
We then take the mean over the $N_{2}$ samples for our estimate of
$D_{\mathrm{KL}}(p(y|\theta_{i},\xi)\|p(y|\xi))$ on line 2 of our algorithm.

###  6.3 REIG as Log-Sum-Exp Stabilization

Step 3 of algorithm 1 shows that the convex dual objective function for the
$I_{\epsilon}$ design criterion manifests as a $\lambda\epsilon$-biased and
$\lambda^{-1}$-weighted log-sum-exp combination of the
$-D_{\mathrm{KL}}(p(y|\theta,\xi)\|p(y|\xi))$ estimates.

From the bounds for log-sum-exp operators we have

$\displaystyle\lambda\epsilon-\min_{i}\\{d_{i}\\}$
$\displaystyle\geq\lambda\epsilon+\lambda\log\frac{1}{N}\sum_{i=1}^{N}\exp(-\lambda^{-1}d_{i})$
$\displaystyle>\lambda\epsilon-\min_{i}\\{d_{i}\\}-\lambda\log N.$

When the ambiguity set radius $\epsilon$ is smaller, the optimal $\lambda$ is
larger, and $M_{\epsilon}$ becomes more like the sample mean; when $\epsilon$
is larger, the optimal $\lambda$ is smaller until eventually the $\lambda\geq
0$ constraint becomes binding. In that case the limit as $\lambda\to 0$ is
achieved and the value is squeezed to become $M_{\epsilon}=-\min_{i}\ d_{i}$.

We interpret these facts in the following way: $I_{\epsilon}$ tends to bias
the samples in the EIG estimate more towards the smaller values of
$D_{\mathrm{KL}}(p(y|\theta_{i},\xi)\|p(y|\xi))$ the larger $\epsilon$ is. A
well-recognized problem in EIG estimation [Foster et al., 2020] is the
presence of samples where $p(y|\xi)$ is under-estimated, leading to
artificially inflated EIG estimates.

We have until this point consider $I_{\epsilon}$ a design criterion in its own
right, but this biasing behavior suggests that the use of algorithm 1 with an
appropriate choice of $\epsilon$ can also be useful as an estimator for the
original EIG criterion whose bias protects the estimate from sampling error.
This may be a viable approach for stabilizing the computed EIG value when
there are insufficient samples to converge the estimate of $p(y|\xi)$ well in
the nested sampling approaches.

Since $\hat{I}^{\text{VNMC}}$ is an upper bound while $\hat{I}^{\text{ACE}}$
and $\hat{I}^{\text{MINE}}$ are lower bounds for $I$, we define an estimator
$\hat{I}_{\epsilon}^{\text{VNMC}}$ that uses VNMC as the estimator for (18)
and the estimators $\hat{I}_{\epsilon,\max}^{\text{ACE}}$ and
$\hat{I}_{\epsilon,\max}^{\text{MINE}}$ that use ACE and MINE, respectively.
The minimization / maximization in each case acts counter to the bias of the
sampling process, which we measure in the following section.

##  7 Experiments

We perform two types of numerical experiments. In the first type, we test to
see how well $I_{\epsilon}$ performs as a “worst-case” EIG estimate, as
discussed in section 5. In the second type, we compare the previously
developed EIG estimators from section 6.2 with a large number of samples to
$I_{\epsilon}$ with a small number of samples to see if the bias against large
summands described in section 6.3 is as effective as additional samples in
stabilizing the computation.

###  7.1 REIG as a Worst-Case Estimate

Figure 2: Preference test: EIG change with prior perturbation and the
performance of $I_{\epsilon}$

This test uses the Preference model test case from [Foster et al., 2019]. The
parameter $\theta$ is a location parameter with a reference prior $p(\theta)$
that is normally distributed, and the experiments are indexed by locations
relative to $\theta$.

To test how well $I_{\epsilon}$ serves at computing worst-case estimate, we
select as second prior $q(\theta)$ such that $I(p)$ and $I(q)$ look visibly
distinct over the range of possible experiments in fig. 2, and then measure
the KL-divergence between $p(\theta)$ and $q(\theta)$ and take this number to
be $\epsilon$ (in this case, $\epsilon=0.2$). We then compute the
$I_{\epsilon}$ criterion for all experiments as well, and plot them in fig. 2.
(In all of these calculations we use a sampling estimator with a large number
of samples so we can be reasonably sure that the values are converged.) What
we see is that $I_{\epsilon}(p,\xi)$ succeeds at being a lower bound for both
$I(p,\xi)$ and $I(q,\xi)$ in this case. We also see that this is not the
result of a uniform or log-uniform scale reduction: $I_{\epsilon}(p,\xi)$ has
more drastically reduced the gain of some experiments than others, meaning
that $I_{\epsilon}(p,\xi)$ determines a different optimal $\xi^{*}$ than
$I(p,\xi)$.

###  7.2 EIG Stabilization via REIG

We now test the effectiveness of $\hat{I}_{\epsilon}(p,\xi)$’s log-sum-exp
stabilization at producing a stabilized estimation of $I(p,\xi)$. We use the
benchmarks designed by [Foster et al., 2019, Kleinegesse and Gutmann, 2020],
especially three experimental designs which have explicit models for the
likelihood distribution: A/B test, Preference, and Pharmacokinetic model.

####  7.2.1 A/B test

An A/B test [Kohavi et al., 2009, Box et al., 1978] can be used to determine
which experiment in a set of two or more results in the largest information
gain. Foster et al. [2019] introduces the group size selection problem between
groups A and B. We have $n$ experiment participants, and we can select $n_{A}$
participants for group A and $n-n_{A}$ participants for group B. Each
participant is represented as two random variables: the first is measured for
participants in group A, and the second for group B. This A/B test models each
participant’s random variables using a Bayesian linear model,
$y=X\theta+\epsilon$; the prior and likelihood distributions are Gaussian
distributions.

Each of the estimators under consideration uses a neural network to generate
samples that musts be trained. The proposal distribution used by the VNMC and
ACE estimators is trained with $A$ and $\Sigma_{p}$ as in [Foster et al.,
2019]: more details about the proposal distribution can be found there.
Likewise, the parameters $\psi$ of the neural network in
$\hat{I}^{\text{MINE}}$ are trained for each $\theta$ with samples of $y$ from
the given distributions.

In our work, we start from a different reference prior distribution, with its
mean taken to be $[4.46,0]$ instead of $[0,0]$. We find that the VNMC
estimator shows considerably more error for some designs with a small number
of samples $M=30$ (fig. 3, top left), and requires more samples to converge
($M\geq 100$). MINE, on the other hand, doesn’t have big difference when we
have enough samples (fig. 3, bottom left).

In contrast to that larger number of samples, we compute
$\hat{I}_{\epsilon}^{\text{VNMC}}$ and $\hat{I}_{\epsilon,\max}^{\text{ACE}}$
estimators using only $M=30$ posterior samples and
$\hat{I}_{\epsilon,\max}^{\text{MINE}}$ using $1000*10$ samples, with an
increasing series of ambiguity set radii $\epsilon$ (fig. 3, right). The
estimate is almost unaffected by the ambiguity set if $\epsilon=0.001$, but
increasing $\epsilon$ to $0.1$ seems to improve the estimates for the designs
that were highly overestimated without negatively affecting the other
experiments that were already properly estimated.

Figure 3: A/B test convergence for EIG estimators (left column): 100*10
samples from $p(\theta,y)$ and 30/100/1000 samples from
$q_{\phi}(\theta|y,\xi)$ (10,000, 30,000, 50,000 samples for the MINE
estimator); $I_{\epsilon}$ estimators (right column): 100 $\theta$ samples
from $p(\theta)$ and 10 y samples from each $\theta_{i}$ and 30 posterior
samples for the marginal likelihood distribution (1000 $\theta$ and $10$ y
samples for the MINE estimator).

####  7.2.2 Preference

Another benchmark that we use to test the effectiveness of
$I_{\epsilon}(p,\xi)$’s log-sum-exp stabilization is Preference test. The
Preference experiment is designed to understand consumer behavior with a
utility function [Samuelson, 1948, Foster et al., 2019]. The experiment
provides the proposal to subjects and checks their preference as an output to
use in the utility function.

We use the experiment from [Foster et al., 2019]. We start from different
reference prior distribution using $-7.35$ as the mean instead of $0$. The KL
divergence from the original prior distribution is $\epsilon=0.2$. VNMC
estimator in our test shows more error with a small number of samples M = 30
(fig. 4, top left), while ACE estimator shows accurate estimation for the EIG
with small number of marginalization samples (fig. 4, middle left). MINE
estimator, on the other hand, estimate the EIG lower than the true EIG (fig.
4, bottom left). But the optimal experimental set-up are same with the true
EIG.

We further estimate $\hat{I}_{\epsilon}^{\text{VNMC}}$ and
$\hat{I}_{\epsilon,\max}^{\text{ACE}}$ estimator using only M = 30 samples and
MINE, with an increasing series of ambiguity set radii $\epsilon$ (fig. 4,
right). By applying ambiguity set 0.01, the $I$ estimation boundary becomes
tight. If we apply ambiguity set $\epsilon=0.1$, the lower bound (ACE) becomes
higher than the true $I$ while the upper bound (VNMC) becomes lower than the
True $I$ (fig. 4, right). The $\hat{I}_{\epsilon,\max}^{\text{ACE}}$ and
$\hat{I}_{\epsilon}^{\text{VNMC}}$ are no longer the corresponding upper bound
and lower bound. $\hat{I}_{\epsilon,\max}^{\text{MINE}}$ with $\epsilon=0.1$
increase the EIG estimation and change the optimal experimental set up (fig.
4, bottom right).

Figure 4: Preference test convergence for EIG estimators (left column): 100*10
samples from $p(\theta,y)$ and 30/100/1000 samples from
$q_{\phi}(\theta|y,\xi)$ (10,000, 30,000, 50,000 samples for the MINE
estimator); $I_{\epsilon}$ estimators (right column): 100 $\theta$ samples
from $p(\theta)$ and 10 y samples from each $\theta_{i}$ and 30 posterior
samples for the marginal likelihood distribution (1000 $\theta$ samples and 10
y samples for MINE estimator).

####  7.2.3 Pharmacokinetic model

The last benchmark that we use to test the effectiveness of
$I_{\epsilon}(p,\xi)$ is Pharmacokinetic study. The Pharamcokinetic study is
designed to understand how the medicine is absorbed, distributed and
eliminated in the subject’s body, which we can understand with the
compartmental model. We can estimate the compartment model’s parameter by
blood sampling and we want to calculate the optimal blood sampling time. [Ryan
et al., 2014, Kleinegesse and Gutmann, 2020]

We start from a different reference prior distribution with its mean taken to
be $[0.1,\log 0.1\log 20]$ whose KL divergence from the original prior
distribution is 0.1. The accuracy of the estimators $\hat{I}^{\text{VNMC}}$
and $\hat{I}^{\text{ACE}}$ in our tests have increases with M (fig. 5, top
left, middle left) while $\hat{I}^{\text{MINE}}$, on the other hand, does
converge to a substantial underestimate (fig. 5, bottom left), yet identifies
the same optimal experiment.

We further compute $\hat{I}_{\epsilon}^{\text{VNMC}}$ and
$\hat{I}_{\epsilon,\max}^{\text{ACE}}$ estimators using only M = 30 samples
and $\hat{I}_{\epsilon,\max}^{\text{MINE}}$, with an increasing series of
ambiguity set radii $\epsilon$ (fig. 5, right). By applying ambiguity set 0.1,
the lower bound (ACE) becomes higher than the true $I$ while the upper bound
(VNMC) becomes lower than the True $I$ (fig. 5, right). The
$\hat{I}_{\epsilon,\max}^{\text{ACE}}$ and $\hat{I}_{\epsilon}^{\text{VNMC}}$
are no longer the corresponding upper bound and lower bound. MINE estimator
with ambiguity set $\epsilon=0.1$, on the other hand, improve the EIG
estimation so that the estimation is close to the true EIG (fig. 5, bottom
right).

Figure 5: Pharmacokinetic model convergence for EIG estimators (left column):
100*10 samples from $p(\theta,y)$ and 30/100/1000 samples from
$q_{\phi}(\theta|y,\xi)$ (10,000, 30,000, 50,000 samples for the MINE
estimator); $I_{\epsilon}$ estimators (right column): 100 $\theta$ samples
from $p(\theta)$ and 10 y samples from each $\theta_{i}$ and 30 posterior
samples for the marginal likelihood distribution (1000 $\theta$ samples and 10
y samples for MINE estimator).

##  8 Discussion

This work presents an introduction and initial numerical experiments testing
the use of a robust modification of expected information gain as a criterion
for Bayesian model selection. The $I_{\epsilon}$ estimator is designed to both
have rigorously defined approximation properties (as the minimization of a
tangent approximation to the EIG over a convex ambiguity set), and to yield a
practical algorithm in practice (algorithm 1, which post-processes the samples
of previously defined estimators by the solution of a 1D convex optimization
problem).

While our initial results are promising, our understanding of how to apply
this method is not complete at this time. Most of the remaining questions
relate to the choice of the radius $\epsilon$ of the ambiguity set over which
the approximated EIG is taken to be robust. A definite a priori estimate for
$\epsilon$ seems unlikely in most cases.

In an approach analogous to Morozov’s discrepancy principle, a logical choice
of $\epsilon$ would be one such that the implied radius of the ambiguity set
is on the same order as the error that the sampling based estimator (such as
VNMC or ACE) has introduced into the problem. The fact that these two
estimators provide (in expectation) an upper and lower bound for the true EIG
of a design suggests that there may be a way to combine the two estimators
into an a posteriori estimate of the proper choice of $\epsilon$. Further
investigation into this topic is the subject of further research.

###### Acknowledgements.

We would like to thank Alexander Shapiro for several stimulating discussions
on ambiguity set.

## References

  * Bayraksan and Love [2015] Güzin Bayraksan and David K Love.  Data-driven stochastic programming using phi-divergences.  In _The Operations Research Revolution_ , pages 1–19. INFORMS, 2015\.
  * Box et al. [1978] George EP Box, William H Hunter, Stuart Hunter, et al.  _Statistics for experimenters_ , volume 664.  John Wiley and sons New York, 1978.
  * DasGupta and Studden [1991] A DasGupta and WJ Studden.  Robust bayesian experimental designs in normal linear models.  _The Annals of Statistics_ , 19(3):1244–1256, 1991.
  * Embretson and Reise [2013] Susan E Embretson and Steven P Reise.  _Item response theory_.  Psychology Press, 2013.
  * Foster et al. [2019] Adam Foster, Martin Jankowiak, Elias Bingham, Paul Horsfall, Yee Whye Teh, Thomas Rainforth, and Noah Goodman.  Variational bayesian optimal experimental design.  In _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.
  * Foster et al. [2020] Adam Foster, Martin Jankowiak, Matthew O’Meara, Yee Whye Teh, and Tom Rainforth.  A unified stochastic gradient approach to designing bayesian-optimal experiments.  In _International Conference on Artificial Intelligence and Statistics_ , pages 2959–2969. PMLR, 2020.
  * Huan [2010] Xun Huan.  _Accelerated Bayesian experimental design for chemical kinetic models_.  PhD thesis, Massachusetts Institute of Technology, 2010.
  * Kleinegesse and Gutmann [2020] Steven Kleinegesse and Michael U Gutmann.  Bayesian experimental design for implicit models by mutual information neural estimation.  In _International Conference on Machine Learning_ , pages 5316–5326. PMLR, 2020.
  * Kohavi et al. [2009] Ron Kohavi, Roger Longbotham, Dan Sommerfield, and Randal M Henne.  Controlled experiments on the web: survey and practical guide.  _Data mining and knowledge discovery_ , 18(1):140–181, 2009.
  * Owhadi et al. [2015] Houman Owhadi, Clint Scovel, and Tim Sullivan.  On the brittleness of bayesian inference.  _siam REVIEW_ , 57(4):566–582, 2015.
  * Ryan et al. [2016a] EG Ryan, CC Drovandi, JM McGree, and AN Pettitt.  Fully bayesian optimal experimental design: A review.  _International Statistical Review_ , 84:128–154, 2016a.
  * Ryan et al. [2014] Elizabeth G Ryan, Christopher C Drovandi, M Helen Thompson, and Anthony N Pettitt.  Towards bayesian experimental design for nonlinear models that require a large number of sampling times.  _Computational Statistics & Data Analysis_, 70:45–60, 2014.
  * Ryan et al. [2016b] Elizabeth G Ryan, Christopher C Drovandi, James M McGree, and Anthony N Pettitt.  A review of modern computational algorithms for bayesian optimal design.  _International Statistical Review_ , 84(1):128–154, 2016b.
  * Ryan [2003] Kenneth J Ryan.  Estimating expected information gains for experimental designs with application to the random fatigue-limit model.  _Journal of Computational and Graphical Statistics_ , 12(3):585–603, 2003.
  * Samuelson [1948] Paul A Samuelson.  Consumption theory in terms of revealed preference.  _Economica_ , 15(60):243–253, 1948.
  * Shapiro [2017] Alexander Shapiro.  Distributionally robust stochastic programming.  _SIAM Journal on Optimization_ , 27(4):2258–2275, 2017.
  * Shapiro et al. [2021] Alexander Shapiro, Enlu Zhou, and Yifan Lin.  Bayesian distributionally robust optimization.  _arXiv preprint arXiv:2112.08625_ , 2021.
  * Stark [2015] Philip B Stark.  Constraints versus priors.  _SIAM/ASA Journal on Uncertainty Quantification_ , 3(1):586–598, 2015.
  * [19] Ruoyong Yang and James O Berger.  _A catalog of noninformative priors_.

##  Appendix A Experiments Details

###  A.1 A/B test

The reference prior and likelihood for the A/B test were as follows:

$\displaystyle\theta\sim\mathcal{N}(\begin{pmatrix}0\\\
0\end{pmatrix},\begin{pmatrix}10^{2}&0\\\
0&1.82^{2}\end{pmatrix}),y|\theta,\xi\sim\mathcal{N}(X_{\xi}\theta,I)$ (25)

###  A.2 Preference

We use the utility function from [Foster et al., 2019].

$\displaystyle\text{let \ }\xi\in\mathrm{R}$ (26)
$\displaystyle\theta\sim\mathcal{N}(-7.35,20^{2})$ (27)
$\displaystyle\eta|\theta,\xi\sim\mathcal{N}(\xi-\theta,1+|\xi|^{2})$ (28)
$\displaystyle y=f(\eta)$ (29) $\displaystyle
f:\mathrm{R}\rightarrow[\epsilon,1-\epsilon]$ (30) $\displaystyle
x\rightarrow\begin{cases}\epsilon&\text{if }x\leq\text{logit}(\epsilon)\\\
1-\epsilon,&\text{if}x\geq\text{logit}(1-\epsilon)\\\
\frac{1}{1-e^{-x}}&\text{otherwise},\end{cases}$ (31)

###  A.3 Pharmacokinetic model

The prior distribution and noise distributions for the pharmacokinetic model
are defined as follows,

$\displaystyle\begin{pmatrix}k_{a}\\\ k_{e}\\\
V\end{pmatrix}\sim\log\mathcal{N}\begin{bmatrix}\begin{pmatrix}0.1\\\ \log
0.1\\\ \log 20\end{pmatrix},\begin{pmatrix}0.05&0&0\\\ 0&0.05&0\\\
0&0&0.05\end{pmatrix}\end{bmatrix}$ (32) $\displaystyle
y_{t}=\frac{400}{V}\frac{k_{a}}{k_{a}-k_{e}}(\exp{-k_{e}t}-\exp{-k_{a}t})(1+\epsilon_{1t})+\epsilon_{2t},$
(33)

where
$\epsilon_{1t}\sim\mathcal{N}(0,0.01),\epsilon_{2t}\sim\mathcal{N}(0,0.1)$,
and $k_{a}>k_{e}$.