https://github.com/nicolelin13/DynamicCalibrationUMAP
https://github.com/nicolelin13/DynamicCalibrationUMAP/Data

# Beyond Static Calibration: The Impact of User Preference Dynamics on
Calibrated Recommendation

Kun Lin  DePaul UniversityChicagoUSA<EMAIL_ADDRESS>,  Masoud Mansoury
Delft University of TechnologyDelftNetherlands<EMAIL_ADDRESS>,  Farzad
Eskandanian  DePaul UniversityChicagoUSA<EMAIL_ADDRESS>,  Milad Sabouri
DePaul UniversityChicagoUSA<EMAIL_ADDRESS>and  Bamshad Mobasher  DePaul
UniversityChicagoUSA<EMAIL_ADDRESS>
(2024)

###### Abstract.

Calibration in recommender systems is an important performance criterion that
ensures consistency between the distribution of user preference categories and
that of recommendations generated by the system. Standard methods for
mitigating miscalibration typically assume that user preference profiles are
static, and they measure calibration relative to the full history of user’s
interactions, including possibly outdated and stale preference categories. We
conjecture that this approach can lead to recommendations that, while
appearing calibrated, in fact, distort users’ true preferences. In this paper,
we conduct a preliminary investigation of recommendation calibration at a more
granular level, taking into account evolving user preferences. By analyzing
differently sized training time windows from the most recent interactions to
the oldest, we identify the most relevant segment of user’s preferences that
optimizes the calibration metric. We perform an exploratory analysis with
datasets from different domains with distinctive user-interaction
characteristics. We demonstrate how the evolving nature of user preferences
affects recommendation calibration, and how this effect is manifested
differently depending on the characteristics of the data in a given domain.
Datasets, codes, and more detailed experimental results are available at:
https://github.com/nicolelin13/DynamicCalibrationUMAP.

Recommender systems, calibration, preference dynamics

††copyright: acmlicensed††journalyear: 2024††doi:
10.1145/3631700.3664869††conference: The 32nd ACM Conference on User Modeling,
Adaptation and Personalization; July 01–04, 2024; Cagliari, Sardinia,
Italy††booktitle: UMAP ’24: ACM Conference on User Modeling, Adaptation and
Personalization, July 01–04, 2024, Cagliari, Sardinia, Italy††isbn:
978-1-4503-XXXX-X/18/06††ccs: Information systems Recommender systems

##  1\. Introduction

Recommender systems learn from past user preferences in order to predict
future user interests and provide users with personalized suggestions.
Traditionally, these systems are evaluated using accuracy metrics like
precision, recall, and normalized discounted cumulative gain (NDCG), measuring
the relevancy of the recommended items delivered to the users. Optimizing
recommender systems for accuracy, however, may result in an imbalance in the
distribution of preference categories in recommendations: the recommendations
may only reflect a user’s dominant preferences in certain categories and
ignore less dominant interests while still being considered accurate. This
type of imbalance, in the long term, can lead to undesirable effects such as
the ”rabbit hole” effect (O’Callaghan et al., 2015; Zhao et al., 2021b) or
filter bubbles (Nguyen et al., 2014).

To address this phenomenon, the notion of Calibration in machine learning has
been extended to the realm of recommender systems (Steck, 2018). The
calibration metric measures the degree to which a user’s prior preference are
reflected in the recommendations generated by the system. For example, in a
hypothetical movie recommender system, a user whose past preferences suggest
70% interest in Action movies and 30% in Drama movies, should ideally see the
same genre ratios in her set of recommended movies. Calibration is typically
measured as the correlation between the distribution of preference categories
in a user’s profile and that of preference categories in generated
recommendations.

To improve calibration (or reduce miscalibration), most existing approaches
assume that users’ preferences distribution is fundamentally static. To
measure calibrations, these approaches typically utilize the full history of
items in users’ profiles to identify preference distributions. This may result
in presumably calibrated recommendations that do not truly reflect users’
current preferences. The problem is that user preferences evolve over time. If
a user’s profile is collected over a long period of time, only certain
segments of that profile may represent the relevant interests of users.

For example, consider a user interacting with a book platform for a period of
10 years. For the first 7 years, the user was intensely interested in books in
the spy-thriller category. But 3 years ago, her interests shifted primarily to
historical novels and then, in the last two years, to science fiction. Using
the standard approach to measuring calibration, a recommender that recommends
items from all three categories based on the user’s full history might be
considered highly calibrated. However, the fact is that the user may not
perceive these recommendations as reflective of her current preferences. On
the other hand, if a recommender only recommends items from the historical and
science fiction categories (i.e., more current user’s preferences), then the
generated recommendations may be considered to be miscalibrated.

Recognizing that a key contributor to miscalibration may be the dynamic nature
of users’ preferences, we hypothesize that selectively utilizing the most
relevant segments from users’ profiles for training the recommendation model
can result in recommendations that better represent users’ current interests.
Our work attempts to shift the focus from predominantly post-processing
methods for re-calibrating recommendations to incorporating calibration
directly into the training phase of the recommendation process, promising a
more fundamental solution to the challenge of maintaining relevance in the
face of evolving user preferences.

In this preliminary work, we conduct an exploratory analysis of calibration in
recommendation systems at a more granular level, taking into account the
dynamics of user preferences. We empirically identify relevant segments of
users’ profiles through a data-driven study of two different datasets with
different user interaction characteristics. Our goal is to show how the
existing approach of measuring calibration that is based on static user
profiles can be misleading and fails to precisely measure the degree to which
recommendations reflect current user interests. To do this, we propose a
preprocessing simulation process that identifies the most representative
segment of the users’ profile by measuring changes in the calibration metric
as user preferences shift over time.

In this simulation process, we first split the users’ profiles into a number
of segments, with the size of segments determined empirically based on the
frequency of user interactions with items. Then, from the most recent segments
to the oldest, we iteratively combine the segments to obtain subsamples of the
data. Finally, on each subsample, we train a recommendation model and measure
recommendation calibration. This process returns the most relevant segment of
a user’s profile that results in the highest calibration with respect to the
current user preference distribution.

Our initial experiments with two datasets from distinctive domains show that
our proposed simulation process enables measuring calibration more precisely
based on users’ changing preferences. We hope that this preliminary work will
lead to more extensive studies of how the dynamics of user preferences in
different domains should be reflected in the design of calibration mechanisms
in recommender systems.

##  2\. Background and Related Work

We denote $\mathcal{U}=\\{u_{1},...,u_{N}\\}$ and
$\mathcal{I}=\\{i_{1},...,i_{M}\\}$ as the set of $N$ users and $M$ items in
the system. We show a user $u$’s profile as $I_{u}$ which contains all
interacted items by $u$. $R_{u}$ is the recommendation list delivered to $u$
which contains all items recommended to $u$. Each item is assigned one or more
categories $\mathcal{C}=\\{c_{1},...,c_{W}\\}$, with $W$ categories in total.

###  2.1. Calibration Metric in Recommendation

The user’s preference distribution over item categories is represented as a
numerical vector of size $W$ that shows the preference ratio of a user toward
each item category. We denote this vector for a user $u$ as
$\mathcal{P}_{u}=\\{p(c|u)|c\in\mathcal{C}\\}$ where $p(c|u)$ is defined as
follows:

(1) $\vspace{-2pt}p(c|u)=\frac{\sum_{i\in I_{u}}{p(c|i)}}{|I_{u}|}$

where $p(c|i)$ is the fraction that category $c$ represents item $i$ (e.g., if
$i$ is assigned $c_{1}$ and $c_{2}$, then $p(c_{1}|i)=0.5$ and
$p(c_{2}|i)=0.5$). Analogously, we define category distribution on the
recommendation list delivered to user $u$, denoted by
$\mathcal{Q}_{u}\\{q(c|u)|c\in\mathcal{C}\\}$, where $q(c|u)$ is defined as
follows:

(2) $\vspace{-2pt}q(c|u)=\frac{\sum_{i\in R_{u}}{p(c|i)}}{|R_{u}|}$

Given the category distribution over the user’s profile and recommendation
list, miscalibration is computed by measuring the distance between these two
distributions. There are various distance measures to consider for this
computation, but following (Steck, 2018), we use Kullback-Leibler (KL)
divergence in this paper:

(3)
$\vspace{-2pt}MC(\mathcal{P},\mathcal{Q})=KL(\mathcal{P}||\mathcal{Q})=\sum_{c\in\mathcal{C}}{p(c|u)\log{\frac{p(c|u)}{q(c|u)}}}$

###  2.2. Methods for Calibrated Recommendation

Steck (Steck, 2018) initially defined calibration emphasizing its role in
achieving personalization that mirrors a user’s past behavior. This concept
has evolved, intersecting significantly with fairness and diversity
(Abdollahpouri et al., 2020; Wang et al., 2023).

Current calibration methods (Sonboli et al., 2020; Eskandanian and Mobasher,
2020; Abdollahpouri et al., 2023), such as those critiqued by (Naghiaei et
al., 2022; Nazari et al., 2022), often separate personalization from
calibration, leading to a potential reduction in recommendation accuracy.
Additionally, these methods typically assume static user preferences,
neglecting the dynamic nature of user interests, as highlighted by Zhao et al.
(Zhao et al., 2021b).

Building upon the insights from previous studies, this research proposes a
novel approach to address miscalibration in recommender systems. Traditional
calibration methods often rely on entire historical data, where older
interactions may no longer reflect the user’s latest preferences, leading to
recommendations that are out of sync with current interests. By selectively
utilizing the most relevant segments of users’ profiles for training the
model, our approach aims to incorporate calibration directly into the training
phase, taking into account the evolving nature of user preferences.

Figure 1. Experiment workflow for dynamic calibration

##  3\. Simulation Process for Dynamic Calibration

In this section, we describe our pre-processing simulation process for
identifying the most representative segments of users’ profiles in training
the recommendation model that would yield more calibrated recommendations.
Figure 1 illustrates this simulation process as follows:

  1. (1)

Given the user-item interaction matrix containing the profile of all users and
their interacted items, we sort each user $u$’s profile $P_{u}$
chronologically in descending order and split $P_{u}$ into $n$ subprofiles as
$\\{P_{u}^{1},P_{u}^{2},...,P_{u}^{n}\\}$ where $P_{u}^{1}$ contains the most
recent $u$’s interactions and $P_{u}^{n}$ contains the oldest interactions.
This process is shown in part (a) of Figure 1. The choice of $n$ (number of
subprofiles) is a domain-specific parameter. The time window size for creating
the subprofiles should be set based on how frequently users interact with
items. For example, in a short video streaming platform where users interact
with items very frequently, a time window of shorter length, like daily, may
be appropriate. However, in a book recommendation platform, a time window of a
longer period, such as one or more years, might better capture users’ evolving
tastes. We discuss and analyze the choice of the time window and $n$ on
different recommendation domains in the methodology section.

  2. (2)

As shown in part (b) of Figure 1, we create samples of the dataset $D$ by
iteratively and chronologically combining the subprofiles of users from
different time windows as follows:

(4) $\vspace{-2pt}D^{l}=\\{P_{u}^{1}\cup P_{u}^{2}\cup...\cup
P_{u}^{l}|\forall u\in\mathcal{U}\\},\;\;\;\;where\;\;l\leq n$

This will create samples of the dataset (i.e., $D^{1},D^{2},...,D^{n}$) as
shown in part (c) of Figure 1.

  3. (3)

We iteratively pick each sample created in the previous step and evaluate the
recommendation performance built on that sample. As shown in part (d) of
Figure 1, we first tune the recommendation model on each sample of the data
for the best model. Then, in part (e) of Figure 1, with the best model
identified in the previous step, we evaluate the recommendation performance
(i.e., miscalibration and accuracy) on each sample of the dataset.

When considering the different segments of the user’s profile, the most recent
segment of the users’ profiles (i.e., $D^{1}$, $D^{2}$, …, or $D^{n}$) that
yields the lowest miscalibration is the one that is most representative of
users’ current preferences.

##  4\. Methodology

###  4.1. Datasets and Preprocessing

We use two datasets from distinct domains. The short video platform is
characterized by rapid user engagement and high interaction frequencies, where
user preferences evolve quickly. For the book domain, users tend to transition
between different book genres more gradually. To better investigate the
calibration dynamics with the shifts in preference patterns, we preprocess
datasets to focus on active users across time windows of different sizes.
Detailed statistics of processed datasets are provided in Table 1

KuaiRec dataset (Gao et al., 2022) is collected on the short video mobile app
Kuaishou111https://www.kuaishou.com/cn from July to September 2020. Each video
is manually tagged with video categories. With watch_ratio as the user’s
preference signal, we refer to its mean (0.9) as the threshold defining the
user’s interest. Positive interactions ($watch\\_ratio\geq 0.9$) are labeled
as ”1” with others as ”0”. Users in this dataset are very active. 86.4% of
users are active for 60 out of the total 63 days, while 99.9% are active for
more than 50 days. In our future work, we will include an analysis of the
larger matrix from KuaiRec with higher sparsity which reflects a more
realistic interaction scenario.

The book recommendation dataset (Wan and McAuley, 2018; Wan et al., 2019) from
GoodReads222https://mengtingwan.github.io/data/goodreads.html includes user-
book interactions. Each book has its assigned genres from all eight book
genres. To reduce dataset size and sparsity from the original dataset we
focused on interactions associated with more engaged users and frequently
accessed items. The dataset was filtered to include only active users, defined
as those engaging in 20 to 50 annual interactions from year 2010 to 2017.
Subsequently, items with fewer than 1,000 interactions were excluded. The
final step filtered in users with a consistent activity level of more than 4
interactions per year.

Table 1. Statistics of the Preprocessed Datasets. Dataset | #Users | #Items | Category type | #Categories | #Interactions | Sparsity
---|---|---|---|---|---|---
KuaiRec | 1,411 | 3,327 | content tag | 31 | 1,799,403 | 61.6%
GoodReads | 865 | 1,662 | book genre | 7 | 104,325 | 92.7%

###  4.2. Tuning Time Windows Representing User Preference Segments

The time window selection over user profile is critical when creating users’
subprofiles in Figure 1 part (a). In this exploratory work, we are interested
in observing changes in calibration when extending the training time window
sizes. Given that goal, we look for the smallest window size, leading to
sufficient interactions for training and representing user preferences. Due to
differences in user interaction frequency and evolving preferences across
domains, time window sizes vary across domains. Users in the book reading
domain have less frequent interactions and more stable preferences than users
in the short video domain. After tuning the time window sizes for the two
datasets (for KuaiRec, 1/7/14 days; for GoodReads, 3/6/12 months), we picked 1
day for KuaiRec and 0.5 years (6 months) for GoodReads.

###  4.3. Recommendation Models

Recommendations are generated by well-tuned Bayesian Personalized Ranking
(BPR) (Rendle et al., 2012)333Our full preliminary experiments include two
algorithms: BPR and ItemKNN. Given the results between the two algorithms are
similar, we focus on BPR in the main paper while making ItemKNN results
available via https://github.com/nicolelin13/DynamicCalibrationUMAP.,
implemented using the library RecBole (Zhao et al., 2021a). With NDCG@10 as
objective, hyperparameters are tuned via an exhaustive search for learning
rate ranging from 0.0001 to 0.01 and embedding size ranging between 64 to 128.
Models for different time periods are independently tuned as shown in Figure 1
part (d). Also, we repeated the experiments six times to ensure the stability
and accountability of the results. The results for analysis are based on the
aggregated results from the repeated experiments.

###  4.4. User Segment Analysis

Given that users’ varied behavioral patterns in the same dataset may be
related to different degrees of calibration across individuals, we extend our
analysis of calibration dynamics to user segments with different
characteristics. We explore two factors for segmenting users: user activity
frequency (number of user interactions) and category-wise entropy in the user
profile. We simply used a percentile-based approach, segmenting users into
three groups based on each factor. A higher entropy indicates that user
preferences are more evenly distributed across categories, while users with
lower entropy have more intense preferences over a few specific categories.
(Lin et al., 2020) showed user profile entropy is an important user profile
factor positively associated with the miscalibration of the received
recommendation.

##  5\. Results and Discussion

###  5.1. Calibration Dynamics Due to Time Window Selections

Fig. 2 depicts the miscalibration dynamics with the time window changing.
Regardless of domains, we observe that miscalibration can vary based on the
time window selections. In KuaiRec, the effect is more pronounced due to the
highly dynamic nature of user preferences, with the optimal calibration
achieved with a window size of 5 days. It is noteworthy that this window size
also corresponds to the optimal NDCG values. In GoodReads, the calibration
peaked at time window 7, indicating the most recent 3.5 years. In this domain,
user preferences tend to be more stable over time, so there is less variance
in calibration values. We note, however, that in this case, recommendation
accuracy increases as larger profile segments are used in training.

((a))

((b))

Figure 2. Box plot of miscalibration distribution by time window

Box plot of miscalibration distribution by time window

###  5.2. Calibration Dynamics Based on User Characteristics

((a))

((b))

Figure 3. Miscalibration by user segments in KuaiRec. The left is segmented by
activity level, and the right is segmented by profile entropy.

Miscalibration by user segments in KuaiRec. The left is segmented by activity
level, and the right is segmented by profile entropy.

((a))

((b))

Figure 4. Miscalibration by user segments in GoodReads. The left is segmented
by activity level, and the right is segmented by profile entropy.

Miscalibration by user segments in GoodReads. The left is segmented by
activity level, and the right is segmented by profile entropy.

For our analysis of the impact of user segment characteristics on calibration
dynamics, we empirically obtain the optimal time window for calibration as
before, but we do so for each user segment instead of the full set of users.
Fig. 3, depicts the changes in miscalibration as window size increases for the
three user segments (with the left panel showing results for user segments
based on interaction frequency and the right panel depicting results for user
segments based on profile entropy). These results suggest that in this domain
there are no significant variances in calibration across different user
segments, especially when considering smaller time windows. This is likely due
to the highly dynamic and unstable nature of user preferences. As before, we
still see an optimal window size of 5 days for all user segments.

Fig. 4 shows a very different picture for the GoodReads domain than that of
KuaiRec. The optimal time window across user segments is not exactly
consistent, but within the overall range that we saw for all users.
Nevertheless, the results indicate that in this domain, it makes sense to
consider different window sizes for distinct user segments to achieve optimal
calibration. However, the more notable observation is that for user segments
with a high activity level, a higher overall calibration is achieved
regardless of window size. We conjecture that this is because users with
higher activity in this domain (users who read books often) tend to have more
stable profiles with little or no change in preference patterns over time.
This is also the case for user segments with low profile entropy. These are
users who tend to have very focused interests in a few genres, and so their
preference patterns tend to be quite stable.

Overall, these results support the conjecture that the degree of calibration
in recommendation depends on the dynamics of user preferences, with more
stable preference profiles resulting in a higher degree of calibration.

##  6\. Conclusion and Future Work

In this exploratory analysis, we have shown that the calibration mechanism in
recommender systems must take into account the dynamic nature of user
preferences. Our simulation-based analysis shows that the degree of
miscalibration decreases as the training time window size changes, focusing on
the most relevant segments of the users’ profiles. Our analysis also shows
that user interaction characteristics associated with different domains lead
to different optimal training window sizes to achieve better calibration.

In future work, we extend this preliminary analysis to include comparison with
post-processing calibration approaches, and also to explore trade-offs between
dynamic calibration and other beyond accuracy metrics such as diversity and
fairness. We will also conduct an analytical investigation on the impact of
domain and user preference patterns on miscalibration dynamics. Currently, we
have only considered two domains with a limited approach to segmenting
usersand provided preliminary conclusions based on these observations. In
future work, we will do a more comprehensive investigation involving
additional data sets with distinct characteristics. We will also perform an
analysis of how different classes of recommendation algorithms, including
those based on session-based, context-aware, and sequence-aware recommendation
models, perform using the dynamic approach to caliberation.

## References

  * (1)
  * Abdollahpouri et al. (2020) Himan Abdollahpouri, Masoud Mansoury, Robin Burke, and Bamshad Mobasher. 2020.  The connection between popularity bias, calibration, and fairness in recommendation. In _Proceedings of the 14th ACM Conference on Recommender Systems_. 726–731.
  * Abdollahpouri et al. (2023) Himan Abdollahpouri, Zahra Nazari, Alex Gain, Clay Gibson, Maria Dimakopoulou, Jesse Anderton, Benjamin Carterette, Mounia Lalmas, and Tony Jebara. 2023.  Calibrated recommendations as a minimum-cost flow problem. In _Proceedings of the Sixteenth ACM International Conference on Web Search and Data Mining_. 571–579.
  * Eskandanian and Mobasher (2020) Farzad Eskandanian and Bamshad Mobasher. 2020.  Using stable matching to optimize the balance between accuracy and diversity in recommendation. In _Proceedings of the 28th ACM Conference on User Modeling, Adaptation and Personalization_. 71–79.
  * Gao et al. (2022) Chongming Gao, Shijun Li, Wenqiang Lei, Jiawei Chen, Biao Li, Peng Jiang, Xiangnan He, Jiaxin Mao, and Tat-Seng Chua. 2022.  KuaiRec: A Fully-Observed Dataset and Insights for Evaluating Recommender Systems. In _Proceedings of the 31st ACM International Conference on Information & Knowledge Management_ (Atlanta, GA, USA) _(CIKM ’22)_. 540–550.  https://doi.org/10.1145/3511808.3557220
  * Lin et al. (2020) Kun Lin, Nasim Sonboli, Bamshad Mobasher, and Robin Burke. 2020.  Calibration in collaborative filtering recommender systems: a user-centered analysis. In _Proceedings of the 31st ACM Conference on Hypertext and Social Media_. 197–206.
  * Naghiaei et al. (2022) Mohammadmehdi Naghiaei, Hossein A Rahmani, Mohammad Aliannejadi, and Nasim Sonboli. 2022.  Towards confidence-aware calibrated recommendation. In _Proceedings of the 31st ACM International Conference on Information & Knowledge Management_. 4344–4348.
  * Nazari et al. (2022) Zahra Nazari, Praveen Chandar, Ghazal Fazelnia, Catherine M Edwards, Benjamin Carterette, and Mounia Lalmas. 2022.  Choice of implicit signal matters: Accounting for user aspirations in podcast recommendations. In _Proceedings of the ACM Web Conference 2022_. 2433–2441.
  * Nguyen et al. (2014) Tien T Nguyen, Pik-Mai Hui, F Maxwell Harper, Loren Terveen, and Joseph A Konstan. 2014.  Exploring the filter bubble: the effect of using recommender systems on content diversity. In _Proceedings of the 23rd international conference on World wide web_. 677–686.
  * O’Callaghan et al. (2015) Derek O’Callaghan, Derek Greene, Maura Conway, Joe Carthy, and Pádraig Cunningham. 2015.  Down the (white) rabbit hole: The extreme right and online recommender systems.  _Social Science Computer Review_ 33, 4 (2015), 459–478.
  * Rendle et al. (2012) Steffen Rendle, Christoph Freudenthaler, Zeno Gantner, and Lars Schmidt-Thieme. 2012.  BPR: Bayesian personalized ranking from implicit feedback.  _arXiv preprint arXiv:1205.2618_ (2012).
  * Sonboli et al. (2020) Nasim Sonboli, Farzad Eskandanian, Robin Burke, Weiwen Liu, and Bamshad Mobasher. 2020.  Opportunistic multi-aspect fairness through personalized re-ranking. In _Proceedings of the 28th ACM Conference on User Modeling, Adaptation and Personalization_. 239–247.
  * Steck (2018) Harald Steck. 2018.  Calibrated recommendations. In _Proceedings of the 12th ACM conference on recommender systems_. 154–162.
  * Wan and McAuley (2018) Mengting Wan and Julian McAuley. 2018.  Item recommendation on monotonic behavior chains. In _Proceedings of the 12th ACM conference on recommender systems_. 86–94.
  * Wan et al. (2019) Mengting Wan, Rishabh Misra, Ndapa Nakashole, and Julian McAuley. 2019.  Fine-grained spoiler detection from large-scale review corpora.  _arXiv preprint arXiv:1905.13416_ (2019).
  * Wang et al. (2023) Chenyang Wang, Yankai Liu, Yuanqing Yu, Weizhi Ma, Min Zhang, Yiqun Liu, Haitao Zeng, Junlan Feng, and Chao Deng. 2023.  Two-sided Calibration for Quality-aware Responsible Recommendation. In _Proceedings of the 17th ACM Conference on Recommender Systems_. 223–233.
  * Zhao et al. (2021a) Wayne Xin Zhao, Shanlei Mu, Yupeng Hou, Zihan Lin, Yushuo Chen, Xingyu Pan, Kaiyuan Li, Yujie Lu, Hui Wang, Changxin Tian, et al. 2021a.  Recbole: Towards a unified, comprehensive and efficient framework for recommendation algorithms. In _proceedings of the 30th acm international conference on information & knowledge management_. 4653–4664.
  * Zhao et al. (2021b) Xing Zhao, Ziwei Zhu, and James Caverlee. 2021b.  Rabbit holes and taste distortion: Distribution-aware recommendation with evolving interests. In _Proceedings of the Web Conference 2021_. 888–899.


# Combining Normalizing Flows and Quasi-Monte Carlo

Charly Andral
<EMAIL_ADDRESS>CEREMADE, CNRS, Université Paris-Dauphine,
Université PSL, 75016 PARIS, FRANCE

###### Abstract

Recent advances in machine learning have led to the development of new methods
for enhancing Monte Carlo methods such as Markov chain Monte Carlo (MCMC) and
importance sampling (IS). One such method is normalizing flows, which use a
neural network to approximate a distribution by evaluating it pointwise.
Normalizing flows have been shown to improve the performance of MCMC and IS.
On the other side, (randomized) quasi-Monte Carlo methods are used to perform
numerical integration. They replace the random sampling of Monte Carlo by a
sequence which cover the hypercube more uniformly, resulting in better
convergence rates for the error that plain Monte Carlo. In this work, we
combine these two methods by using quasi-Monte Carlo to sample the initial
distribution that is transported by the flow. We demonstrate through numerical
experiments that this combination can lead to an estimator with significantly
lower variance than if the flow was sampled with a classic Monte Carlo.

##  1 Introduction

Sampling from a distribution known only by its (unnormalized) density is a
problem that has arisen in many different fields such as physics, Bayesian
statistics and machine learning. In particular it if often used to estimate an
expectation of a function with respect to this distribution, which is known as
Monte Carlo integration [29]. Monte Carlo has been a very important topic in
statistics in recent decades [42], especially with the development of Markov
Chain Monte Carlo (MCMC), as the Metropolis-Rosenbluth-Teller-Hastings (MRTH)
algorithm [28, 16] or importance sampling (IS) [21, 12].

In the recent years, generative methods coming from the machine learning
community has been used to enhance existing Monte Carlo methods [19, 41]. In
particular, normalizing flows [46, 45, 8] offer a general framework to
approximate the distribution of interest by a measure parametrized by a neural
network, while being able to sample from it and compute its density. The
samples from the flow can then be used as the instrumental distribution of an
importance sampler [31, 32] or as a proposal for a MRTH kernel [1, 11].
Recently, a comparison of the different ways to use flow to sample was done
[14].

On the other hand, (randomized) quasi-Monte Carlo methods [35] have been for
decades to perform numerical integration by replacing the random sample of
Monte Carlo by a sequence that covers the hypercube $[0,1]^{d}$ more
uniformly, allowing to obtain better convergence rates than the $O(n^{-1/2})$
of plain Monte Carlo, typically at $O(n^{-1}(\log(n)^{d}))$ for QMC. Quasi-
Monte Carlo can be combined with MCMC [39, 5, 7, 26] or with importance
sampling [6, 22].

As far as we are aware of, this combination of normalizing flows with quasi-
Monte Carlo has only been tried in [47]. The scope of the paper was different,
as the authors were interested in QMC to minimize the error of the stochastic
gradient for the training of the flow, even though they mentioned that the QMC
sample could be used to reduce the variance when estimating an expectation. In
particular, the combination with IS and MCMC was not studied.

In this work, we explore the combination of using normalizing flows as an
approximation of the target distribution and using quasi-Monte Carlo to sample
the initial distribution of the flow. We will show that in same cases, this
combination lead to a estimator with a significantly lower variance than if
the flow was sampled from i.i.d variables.

##  2 Normalizing Flows

###  2.1 Using an instrumental distribution for Monte Carlo

We place ourselves in the classic Monte Carlo setting, where, for a given
distribution $\pi$ on $\mathbb{R}^{d}$, we want to estimate the expectation of
a test function $f$ with respect to $\pi$, defined as

$\pi(f):=\mathbb{E}_{\pi}\left[f(X)\right]=\int_{\mathbb{R}^{d}}f(x)\pi(dx)$

in the case where a direct sample from $\pi$ is not available and the integral
is intractable.

##### Importance Sampling

One solution, if we can sample from an instrumental distribution $\rho$ close
enough to $\pi$, is to use the importance sampling estimator by noticing that:

$\pi(f)=\int_{\mathbb{R}^{d}}f(x)\frac{\pi(x)}{\rho(x)}\rho(dx)=\mathbb{E}_{\rho}\left[f(X)\frac{\pi(X)}{\rho(X)}\right].$

Therefore, if we can sample from $\rho$, we can estimate $\pi(f)$ by:

$\pi(f)\approx\hat{I}_{\text{IS}}(f):=\frac{1}{N}\sum_{i=1}^{N}f(X_{i})w(X_{i}),$
(1)

for $X_{i}\sim\rho$ and $w(x)=\frac{\pi(x)}{\rho(x)}$, which is called the
importance weight.

If the densities $\pi$ and/or $\rho$ are known only up to a normalizing
constant, the weight $w(x)$ is known up to a constant. We denote the
unnormalized importance weight by $\tilde{w}$ := $c\cdot w$ for an unknown
$c\in\mathbb{R}$. In that case, the estimator $\hat{I}_{\text{IS}}(f)$ cannot
be used as is and is replace by the self-normalized importance sampling
estimator:

$\hat{I}_{\text{SNIS}}=\frac{\sum_{i=1}^{n}\tilde{w}(X_{i})f(X_{i})}{\sum_{i=1}^{n}\tilde{w}(X_{i})}.$
(2)

While $\hat{I}_{\text{IS}}$ is unbiased (i.e.
$\mathbb{E}[\hat{I}_{\text{IS}}]=\pi(f)$), $\hat{I}_{\text{SNIS}}$ is biased
due to the ratio of the two estimates.

##### Independent MRTH algorithm

Another solution is to use a Metropolis-Rosenbluth-Teller-Hastings (MRTH)
kernel using $Q(x,dy)=\rho(dy)$ as the proposal kernel, which we call the
independent MRTH (iMRTH) kernel. This produces a Markov chain that admits
$\pi$ as its invariant distribution. More precisely, at a given time $t$ and
point $X_{t}$, the new state $X_{t+1}$ is defined by

$X_{t+1}=\begin{cases}Y_{t+1}&\text{with probability }\alpha(X_{t},Y_{t+1})\\\
X_{t}&\text{otherwise }\end{cases},$

where $Y_{t+1}$ is drawn from $\rho$ and
$\alpha(x,y):=\min\left(1,\frac{\pi(y)\rho(x)}{\rho(y)\pi(x)}\right)$ is
called the acceptance rate.

##### Importance Markov Chain

Recently [2] proposed a new method called Importance Markov Chain (IMC) that
is in between importance sampling and MCMC. In its independent version (iIMC),
it corresponds to drawing a sample $X_{i}\sim\rho$ and $N_{i}\in\mathbb{N}\sim
R(X_{i},.)$ for $R$ a kernel such that
$\mathbb{E}\left[N_{i}\middle|X_{i}\right]=\kappa\cdot w(X_{i})$. The
parameter $\kappa$ is a fixed constant and $w$ the importance weight function.
The final chain is obtained by concatenating over $i\in 1:n$ the $X_{i}$
replicated $N_{i}$ times.

###  2.2 Normalizing Flows

Normalizing flows are a type of model that can transform a simple distribution
into a more complex one using a neural network trained to target a
distribution through pointwise evaluation.

Formally, a normalizing flow is an invertible and differentiable map $T$ from
$\mathbb{R}^{d}$ to $\mathbb{R}^{d}$ . Given a reference probability measure
$\mu$ (often a standard Gaussian) we can define a new probability measure
$\nu$ as the pushforward of $\mu$ by $T$, i.e.
$\nu=T_{\sharp}\mu=\mu(T^{-1}(.))$. The flow is trained by minimizing the
Kullback-Leibler divergence between $\nu$ and the target distribution $\pi$. A
sample from the flow can be generated by sampling $X\sim\mu$ and then taking
$T(X)\sim\nu$.

The density of the pushforward measure $\nu$ can be computed using the change
of variable formula:

$\nu(x)=\mu(T^{-1}(x))\left|\det(J_{T^{-1}}(x))\right|,$

where $J_{T^{-1}}$ is the Jacobian of $T^{-1}$, the inverse map of $T$.

Normalizing flows are designed to satisfy several useful properties: the
inverse map $T^{-1}$ and the determinant of its Jacobian should be easily
evaluated to allow for the computation of the density of $\nu$; $T$ should be
easily computable to efficiently generate samples from $\nu$; and the class of
flows should be flexible enough to approximate any distribution. For a more
detailed introduction of normalizing flows and the different models, we refer
to [23, 40].

In our numerical experiments, we will use two different types of flows: the
real-valued non-volume preserving (RealNVP) transforms of [9] and the rational
quadratic neural spline flows (RQ-NSF, in their coupling form) of [10], with
the latter being more flexible than the former.

##### Coupling flows

Coupling flows, introduced by [8], are used in RealNVP and RQ-NSF. Each layer
of the flow is defined by a function $f$ that takes an input vector
$x=(x_{i})\in\mathbb{R}^{d}$, a $d_{1}<d$, and defines the output $y=f(x)$
elementwise as follows:

$f(x_{i})=\begin{cases}x_{i}&\text{if }i\leq d_{1}\\\
g(x_{i},\theta(x_{1:d_{1}}))&\text{otherwise}\end{cases}.$ (3)

Here, $g$ is a certain function with a parameter $\theta$ that depends only on
$x_{1:d_{1}}:=(x_{1},\dots,x_{d_{1}})$. For RealNVP, $g$ is an affine
function; for RQ-NSF, $g$ is a rational quadratic spline [13]. The parameter
$\theta$ is called the conditioner and is usually a neural network (multilayer
perceptron) with $x_{1:d_{1}}$ as input. In practice, $d_{1}=d/2$ and the
function uses a mask to alternate the dimensions that are kept unchanged and
the ones that are transformed.

If $g$ is invertible, $f$ is invertible too, and the inverse of $f$ is given
by:

$f^{-1}(y_{i})=\begin{cases}y_{i}&\text{if }i\leq d_{1}\\\
g^{-1}(y_{i},\theta(y_{1:d_{1}}))&\text{otherwise}\end{cases}.$

By construction, the Jacobian of $f$ (and $f^{-1}$) is triangular, and
therefore its determinant is the product of the diagonal elements.

For a more general definition of coupling flows where $f$ is not applied
elementwise and for the defintion of autoregressive flows (not used in this
article), we refer to [23].

##### Training of the flow

In most of the cases (e.g. for RealNVP and RQ-NSF), the flow is parametrized
by a neural network. Let us write $T_{\theta}$ for $\theta\in\mathbb{R}^{d}$
and $\nu_{\theta}:=(T_{\theta})_{\sharp}\mu$ to explicit the dependence of the
flow on the parameter $\theta$. The goal of the training in to find $\theta$
such that $\nu_{\theta}$ is close to $\pi$.

This proximity is usually measured by the (forward) Kullback-Leibler (KL)
divergence between the target distribution $\pi$ and the distribution of the
flow $\nu_{\theta}$:

$D_{\mathrm{KL}}\left(\mu\,\middle\|\,\nu_{\theta}\right):=\int_{\mathbb{R}^{d}}\log\left(\frac{\mu(x)}{\nu_{\theta}(x)}\right)\mu(x)dx.$
(4)

In some cases, the reverse KL divergence
$D_{\mathrm{KL}}\left(\mu\,\middle\|\,\nu_{\theta}\right)$ is minimized
instead. For more details on the differences between the forward and reverse
KL divergences see [40].

This paper does not propose a new training for the flow. Therefore we will use
different trainings from other works in our numerical experiments.

##  3 (Randomized) Quasi-Monte Carlo

###  3.1 Quasi-Monte Carlo

Quasi-Monte Carlo (QMC) methods are used to solve the problem of integrating a
function over the unit hypercube $[0,1]^{d}$. Unlike (plain) Monte Carlo
methods, QMC methods use a deterministic sequence of points that is designed
to cover more uniformly the space, rather than random points. For an
introduction to QMC, see [35]. The uniformity of the sequence is measured by
the star-discrepancy, defined fora set of points
$x_{1},\dots,x_{n}\in[0,1]^{d}$ by:

$D^{\star}_{n}(x_{1},\dots,x_{n})=\sup_{\mathbf{a}\in[0,1]^{d}}\left|\frac{1}{n}\sum_{i=1}^{n}\mathbbm{1}_{x_{i}\in[0,\mathbf{a})}-\prod_{i=1}^{d}a_{i}\right|,$

where $\mathbf{a}=(a_{1},\dots,a_{d})$ and
$[0,\mathbf{a})=\prod_{i=1}^{d}[0,a_{i})\subseteq[0,1]^{d}$.

The error of the QMC estimator built with $x_{1},\dots,x_{n}$ is controlled by
the Koksma-Hlawka inequality [24, 18]:

$\left|\frac{1}{n}\sum_{i=1}^{n}f(x_{i})-\int_{[0,1]^{d}}f(x)dx\right|\leq
D^{\star}_{n}(x_{1},\dots,x_{n})V_{\text{HK}}(f),$

where $V_{\text{HK}}(f)$ is the total variation of $f$ in the sense of Hardy
and Krause. For $d=1$, it coincides with the classic total variation. For the
general definition, see [35].

A QMC sequence is constructed such that
$D^{\star}_{n}(x_{1},\dots,x_{n})=O(n^{-1}(\log n)^{d})$. Therefore, if
$V_{\text{HK}}(f)<\infty$, the QMC estimator has a convergence rate of
$O(n^{-1}(\log n)^{d})$, which is better than the $O(n^{-1/2})$ of plain Monte
Carlo. There are two types of QMC methods: lattice rules and sequences.

##### Lattice rules

Lattice rules are the first type of QMC methods. They involve taking a
generating vector $z\in\mathbb{N}^{d}$ and a fixed $n$, and defining a rank-1
lattice by:

$x_{i}=\frac{(i-1)z}{n}\mod 1,$

for $i\in 1:n$, where the modulo is taken componentwise. The generating vector
$z$ must be carefully chosen to have a low discrepancy. To be efficient,
exactly $n$ points must be used for the estimate, and more points are needed,
a new lattice must be computed. Lattice rules are particularly used in
numerical analysis, but not much in statistics. We will in this paper focus on
the second type of QMC methods, the sequences and digital nets.

##### QMC sequences

The second type of QMC methods are sequences, such as the Halton sequence or
the Sobol’ sequence. [15, 44]. The latter is part of what is called digital
nets. Unlike lattice rules, sequences are infinite, and more points can be
added to the sequence without having to recompute the whole sequence. However,
due to the construction of the Sobol’ sequence, to keep the balance of the
sequence, it is highly recommended to use a number of points $n$ that is a
power of 2, and not to drop the first point [36]

###  3.2 Randomized Quasi-Monte Carlo

As a QMC sequence is deterministic, it does not allow for the estimation of
the error of integration. Therefore, a variant of quasi-Monte Carlo called
randomized quasi-Monte Carlo (RQMC) was developed. The sequence
$x_{1},\dots,x_{n}$ is no longer deterministic, but still has an almost surely
$O(n^{-1}(\log n)^{d})$ star-discrepancy while having every marginal
distribution uniform on $[0,1]^{d}$. The latter property also implies that the
estimator is unbiased, while the plain QMC estimator is usually biased. For
instance, one can take a lattice rule $x_{1},\dots,x_{n}$and apply a common
shift $U\sim\mathcal{U}([0,1]^{d})$ to all the points: $\tilde{x}_{i}=x_{i}+U$
[4].

For sequences like the Halton or Sobol’ sequences, the randomization is
different. Instead of a global shift, the randomization is done by scrambling
the sequence. In a nutshell, scrambling can be describe as a permutation of
the digits (in a certain base) of the points. The permutations are the same
for all the points. Specific scrambling was developpd for the Halton sequence
[33] and the Sobol’ sequence [37, 27].

In the case of scrambled Sobol’ sequence, the randomization can improve the
error from $O(n^{-1}(\log n)^{d-1})$ to $O(n^{-3/2}(\log n)^{(d-1)/2})$ [38].
Running RQMC instead of QMC allows for the estimation of the error by
repeating the experiment and even computing confidence intervals [25].

##  4 Quasi-Monte Carlo for Normalizing Flows

###  4.1 From $[0,1]^{d}$ to $\mathbb{R}^{d}$

In some cases, an integral on $\mathbb{R}^{d}$ can be transformed via a map as
an integral on $[0,1]^{d}$. For instance, a multidimensional Gaussian random
variable can be written as a transform $S:[0,1]^{d}\to\mathbb{R}^{d}$ of a
uniform random variable on the hypercube. This can be done by taking the
inverse transform or the Box-Muller transform [3]. This allows for the use of
QMC.

The inverse transform is defined as follows: if $U_{1}\sim\mathcal{U}([0,1])$
and $F$ is the cumulative distribution function of a standard Gaussian, then
$F^{-1}(U_{1})\sim\mathcal{N}(0,1)$. A multidimensional standard Gaussian
random variable can be obtained by applying this transform componentwise to
$U\in\mathcal{U}([0,1]^{d})$.

The Box-Muller transform is defined as follows: if
$U_{1},U_{2}\sim\mathcal{U}([0,1])$, then $X_{1}=\sqrt{-2\log(U_{1})}\cos(2\pi
U_{2})$ and $X_{2}=\sqrt{-2\log(U_{1})}\sin(2\pi U_{2})$ are two independent
standard Gaussian variables. This can be extended in higher dimension by
taking $U\in\mathcal{U}([0,1]^{d})$ and applying the transform by blocks of
size 2.

The (R)QMC sample $x_{1},\dots,x_{n}$ is now transformed to the Gaussian RQMC
sample $S(x_{1}),\dots,S(x_{n})$. It is worth noting that in that case, as
this transform is not bounded, $V_{\text{HK}}(S)=\infty$ and the Koksma-Hlawka
inequality is useless. However, RQMC can still show better convergence rate
than plain MC (and QMC) [34].

###  4.2 Importance sampling with NF and QMC

First, let us rewrite $\pi(f)$ as an integral over the hypercube $[0,1]^{d}$.
From now on, the instrumental distribution for the IS is set to be the measure
of the flow $\rho$. Recall that if $X \sim\mathcal{N}(0,I_{d})$, then
$T(X)\sim\rho$. Denote S as the mapping from $[0,1]^{d}$ to $\mathbb{R}^{d}$
such that if $U\sim\mathcal{U}([0,1]^{d})$, then $S(U)\sim\mathcal{N}(0,1)$.
Then, by composing the two maps, we have that $(T\circ S)(U)\sim\rho$. We can
finally write $\pi(f)$ as an integral on the hypercube:

$\pi(f)=\int_{[0,1]^{d}}\left(f\cdot w\right)\left((T\circ S)(u)\right)du$

where $w(x)=\pi(x)/\rho(x)$ the importance weight.

Taking $x_{1}, \dots,x_{n}$ from a RQMC sequence, we get the following
estimator:

$\hat{I}^{\text{RQMC}}_{\text{IS}}(f)=\frac{1}{n}\sum_{i=1}^{n}\left(f\cdot
w\right)\left((T\circ S)(x_{i})\right).$ (5)

If the densities are not normalized, the (self-normalized) estimator is :

$\hat{I}^{\text{RQMC}}_{\text{SNIS}}(f)=\frac{\sum_{i=1}^{n}\left(f\cdot\tilde{w}\right)\left((T\circ
S)(x_{i})\right)}{\sum_{i=1}^{n}\tilde{w}((T\circ S)(x_{i}))}$ (6)

As marginally $x_{i}\sim\mathcal{U}(0,1)$ for all $i\in 1:n$,
$\hat{I}^{\text{RQMC}}_{\text{IS}}(f)$ is an unbiased estimator of $\pi(f)$
and $\hat{I}^{\text{RQMC}}_{\text{SNIS}}(f)$ is biased.

Using the exact same notation, but using a classic Monte Carlo sample
$x_{1},\dots,x_{n}$ i.i.d. and uniformly distributed on $[0,1]^{d}$, we can
define the two associated Monte Carlo estimators
$\hat{I}^{\text{MC}}_{\text{IS}}(f)$ and
$\hat{I}^{\text{MC}}_{\text{SNIS}}(f)$. They coincide with the estimator from
(5) and (6) if we do the change of variables.

The case of importance sampling using QMC and Gaussian instrumental
distribution was recently studied in [17] and it can achieve a
$O(n^{-1+\varepsilon})$ converge rate in some cases. A study of the star-
discrepancy in the specific case of importance sampling was conducted in [6].

##  5 Numerical experiments

For the different experiments, we will use two flows from the literature and
show how using quasi-Monte Carlo instead of Monte Carlo can lead to a use
variance reduction. In all the experiments, we will use the same number of
points for the RQMC and the MC estimator. The RQMC sample are computed using
the SciPy package in Python [43]. The implementation for the Sobol’ sequence
use the work of [20] and the scrambling of [27]. The implementation for the
Halton sequence follows [33]. No pure QMC was used; every QMC mentioned in the
figures is RQMC. The code for the experiments is available at
https://github.com/charlyandral/QMC_norm_flows.

###  5.1 Mixture of Gaussians in 2D

####  5.1.1 Setting

For this first experiment, we take the mixture of Gaussians example of [30],
following the code available at https://github.com/lollcat/fab-torch. The
target density $\pi$ is mixture of 40 Gaussians on $\mathbb{R}^{2}$. The
training is done using an $\alpha$-divergence for $\alpha=2$:

$D_{2}(\pi\|\nu_{\theta}):=\int\frac{\pi(x)^{2}}{\nu_{\theta}(x)}dx=\mathbb{E}_{g}\left[\frac{\pi(X)^{2}}{\nu_{\theta}(X)g(X)}\right]$

where $g$ is an importance function. They took $g$ to be an approximation of
$\pi^{2}/\nu_{\theta}$ using annealed importance sampling. The flow is a
realNVP from [9] with 15 layers.

In this example, as the normalizing constant of the target is known, the (non-
self-normalized) importance sampling estimator can be use as defined by (5).

For a same flow, we compare the value of the importance sampling estimator
using a classic Monte Carlo sample and a scrambled Sobol’ sequence with an
inverse map. A total of $2^{16}$ points are taken for every estimate. This is
repeated 200 times to get estimate of the variance and the boxplot. Seven
different test functions are evaluated : $\phi_{1}:x\mapsto x_{1}$,
$\phi_{2}:x\mapsto x_{2}$, $\phi_{3}:x \mapsto x_{1}^{2}$, $\phi_{4}:x\mapsto
x_{2}^{6}$, $\phi_{5}:x:x_{1}x_{2}$,
$\phi_{6}:x\mapsto\sin(x_{1})\cos(-x_{2}/10)$,
$\phi_{7}=\mathbf{1}_{x_{1}>30}$.

####  5.1.2 Results

The results can be visualized in the boxplot of Figure 1, which shows the
different values of $\hat{I}^{\text{RQMC}}_{\text{IS}}(\phi_{i})$ and
$\hat{I}^{\text{MC}}_{\text{IS}}(\phi_{i})$. For the sake of readability, the
values of the estimator are normalized for every function, i.e. the mean of
the estimates is first subtracted and then the result is divided by the
standard deviation. As the this linear scaling is the same for the two
estimators (for a fixed $\phi_{i}$), the boxplots are still comparable.

The ratio of the empirical standard deviations of the two estimators

$r_{i}=\frac{\hat{\sigma}\left(\hat{I}^{\text{MC}}_{\text{IS}}(\phi_{i})\right)}{\hat{\sigma}\left(\hat{I}^{\text{RQMC}}_{\text{IS}}(\phi_{i})\right)}$
(7)

is given in Table 1.

Figure 1: Comparison of the importance sampling estimator using an i.i.d. sample and a Sobol’ sequence for different test functions for the Gaussian mixture. Table 1: Ratio of the standard deviations between RQMC and MC estimates for several test functions for the gaussian mixture. $\phi_{1}$ | $\phi_{2}$ | $\phi_{3}$ | $\phi_{4}$ | $\phi_{5}$ | $\phi_{6}$ | $\phi_{7}$
---|---|---|---|---|---|---
2.829 | 2.285 | 2.182 | 2.787 | 2.757 | 2.465 | 2.374

The results for the seven different test functions are quite similar. They
show a significant variance reduction when using a scrambled Sobol’ sequence
instead of a classic Monte Carlo sample. The ratio of the standard deviations
$r_{i}$ is quite consistent between the different functions, going from
$2.182$ to $2.829$.

###  5.2 Dualmoon target

####  5.2.1 Setting

For the second experiment, we will follow the methodology of [11] to train the
flow. This is implemented in the package FlowMC in Python, [48]. The target we
used is called the dualmoon in this package. In this case, the authors choose
the forward KL divergence (4) as the training objective. As it requires
computing an expectation with respect to the target, direct simulation in not
possible, and the optimization of the flow is done concurrently with an
adaptive MCMC.

The target distribution $\pi$ is defined, for a dimension $d$ that can change,
by:

$\displaystyle\pi(x)\propto$
$\displaystyle\exp\bigg{\\{}-\frac{1}{2}\left(\frac{\left\|x\right\|-2}{0.1}\right)^{2}+$
$\displaystyle\sum_{i=1}^{d}\log\left(e^{-\frac{1}{2}(\frac{x_{i}+3}{0.6})^{2}}+e^{-\frac{1}{2}(\frac{x_{i}-3}{0.6})^{2}}\right)\bigg{\\}}.$

The density is known up to a normalizing constant, and the self-normalized
importance sampling estimator (6) is therefore used for the rest of this
section. The target presents multimodality as the number of modes is $2^{d}$.

The flows are the rational quadratic neural spline flows (RQ-NSF) from [10],
in their masked coupling form. The number of layers and bins will change in
the different illustrations. Except of section 5.2.4, the RQMC sequence is
always a scrambled Sobol’ sequence with a inverse transform.

####  5.2.2 Effects of the dimension

In this first experiment, we study the effect of the dimension on the variance
of the estimator. We will take $d=2,3,4,5,6,7,8,9,10$ and for each dimension,
we train 5 flows with 6 layers and 8 bins and a $(32,32)$ perceptron for the
hidden layers. For each dimension, we take 50 samples of $2^{16}$ points from
scrambled Sobol’ sequences and classic Monte Carlo samples. Two test functions
are evaluated: $\phi_{1}:x\mapsto x_{1}$ and $\phi_{2}:x\mapsto\sin(10\cdot
x_{1})\cdot x_{1}^{4}$. Note that this time, as the target $\pi$ is symmetric
and the two test functions are odd, the true value to estimate is
$\pi(\phi_{1})=\pi(\phi_{2})=0$. This allows computing the ratio if the
absolute errors between the two estimators. Table 2 shows the ratio of the
standard deviations between the two estimators and for the two test functions,
Table 3 shows the ratio of the absolute errors.

dim | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---|---
$\phi_{1}$ | 15.06 | 4.60 | 1.88 | 1.99 | 1.73 | 1.68 | 1.48 | 1.59 | 1.37
$\phi_{2}$ | 39.03 | 6.31 | 2.71 | 2.49 | 2.14 | 1.85 | 1.85 | 1.69 | 1.71
Table 2: Effect of the dimension on the ratio of the standard deviations between RQMC and MC estimates for the dualmoon target. dim | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---|---
$\phi_{1}$ | 20.04 | 4.80 | 2.02 | 2.00 | 1.68 | 1.63 | 1.51 | 1.59 | 1.39
$\phi_{2}$ | 48.31 | 6.42 | 2.79 | 2.61 | 2.21 | 1.88 | 1.92 | 1.69 | 1.72
Table 3: Effect of the dimension on the ratio of the absolute errors between
RQMC and MC estimates for the dualmoon target.

The results in Tables 2 and 3 show very similar results: the efficiency of
RQMC over MC is decreasing with the dimension and also depends on the test
function. The test function $\phi_{1}$ is more simple than $\phi_{2}$, and
thus the variance reduction is more significant. In dimension 2, the
improvement is substantial, with a ratio of standard deviations up to $39.03$
and up to $48.31$ for the ratio of errors. For higher dimensions, the two
types of ratios shows very similar results, less important than in dimension 2
but still significant. As RQMC has better convergence rate that MC, the ratio
is expected to increase with the number of samples. This will be detailed in
section 5.2.4.

This reduction of efficiency with the dimension is not surprising as the
quality of the QMC estimate highly depends on the regularity of the function
to integrate, and the flow becomes more complex with the dimension. In
particular, RQMC is sensitive to effective dimension [35, chapter 17.2].
Intuitively, the effective dimension is a way to measure the complexity of the
function to integrate. The effective dimension if low if the function can be
written a sum of functions depending on a small number of inputs. An RQMC
integration will perform better compared to MC if the effective dimension if
small, even if the dimension of the space $d$ is large. If the two test
functions $\phi_{i}$ seems to have an effective dimension of one, the true
function to integrate with RQMC is $\left(\phi_{i}\cdot w\right)\left((T\circ
S)\right)$. The flow $T$ defined by a coupling of equation 3 is designed to
make the entries interact. Therefore, the effective dimension of the true
function to integrate is higher than the effective dimension of $\phi_{i}$.

####  5.2.3 Effect of the complexity of the flow

We compare the effect of the complexity of the flow on the estimation. For a
fixed dimension $d$, we train five flows in each of the setting form the
parameters, written as (number of layers, hidden size of the conditioner,
number of bins): (2,[8,8],4), (4,[16,16],8), (6,[32,32],8) and
(10,[32,32],12). All flows have the same number of training loops. For each
flow, we take 50 samples of $2^{16}$ points from scrambled Sobol’ sequences
and classic Monte Carlo samples. We evaluate two test functions:
$\phi_{1}:x\mapsto x_{1}$ and $\phi_{2}:x\mapsto\sin(10\cdot x_{1})\cdot
x_{1}^{4}$. The results are shown in Figure 2.

Consistent with the previous section, the variance reduction is more
significant for the simple test function $\phi_{1}$ than for $\phi_{2}$.
Although the simplest flow has a slightly higher average error, the other
three flows have similar results. The most performing is the (6,[32,32],8)
flow. The results obtained here are consistent with the average training loss
of the four flows, respectively 1.43, 0.50, 0.42, 0,43. The ratio of RQMC
error over MC is similar is the four cases.

Figure 2: Impact of the complexity of the flow on absolute error between RQMC
and MC estimates for the dualmoon target in dimension 2

####  5.2.4 Effect of the choice of the RQMC sequence

In this section, we compare Sobol’ and Halton sequences, as well as two
mappings from $\mathcal{U}((0,1)^{d})$ to $\mathcal{N}(0,I_{d})$, the inverse
and the Box-Muller transform described in section 4.1. This gives a total of
four methods. The dimension is set to $d=2$ as it was showed in the previous
section that RQMC is more efficient in low dimensions. This allows for a
better analysis of the choice of the sequence. For visibility, we only show
the results for the two test functions $\phi_{1}:x\mapsto x_{1}$ and
$\phi_{2}:x\mapsto\sin(10\cdot x_{1})\cdot x_{1}^{4}$. For each test function,
we take 100 samples of $2^{17}$ points from the different sequences/transforms
for a common flow. The results are shown in Figure 3.

Figure 3: Convergence of the error of the importance sampling estimator for
the dualmoon target in dimension 2 in log-log scale. True means inverse
transform, false means Box-Muller. Dotted line corresponds to the convergence
rate of $O(n^{-1/2})$ of plain Monte Carlo.

The Figure 3 plots, in log-log scale, the mean of the absolute value of the
error as a function of the sample size. We see that the two graphs are
similar, up to a a vertical shift, corresponding to a multiplicative constant
in the converging rate. The error for MC follow a $O(n^{-1/2})$ convergence
rate, as expected. For RQMC, the four methods give a similar error, with a
convergence rate that accelerate with the sample size (the slope of the line
is increasing). The best performing method is the Sobol’ sequence with a Box-
Muller transform. Overall, Sobol’ sequences perform better than Halton
sequences, and the Box-Muller transform is better than the inverse transform.
Worth noticing that for the Sobol’ sequence, the error is significantly lower
when the sample size is a power of 2. This is coherent with other results on
Sobol’ sequences [36].

####  5.2.5 Results for Markov chains methods

For the last numerical experiment we try two Markov-chain-based methods : the
independent MRTH (iMRTH) and the independent IMC (iIMC) defined in section
2.1. Both methods use i.i.d. proposals from the instrumental distribution
$\rho$ to construct a Markov chain that admits $\pi$ as invariant measure. As
for importance sampling, we replace the $i.i.d.$ sample from $\rho$ for an
RQMC sample using different sequences. The same data are used as in the
previous section, with the same flow. All the analysis is done for the mean of
the first component. The results are shown in Figure 4.

Figure 4: Results for the iMRTH and iIMC methods for the dualmoon target in
dimension 2. (Left) the effective sample size for the different sequences.
(Right) the absolute error of the estimate for the different sequences, as
well as the error of the importance sampling estimator.

The effective sample size is a common tool for Markov chain analysis and is
computed using the autocorrelations the chain [42]. The higher the value, the
better the chain. The plot on the left of Figure 4 shows the effective sample
size for the different sequences. For the iIMC, it shows a great change of
efficiency between the different sequences and transforms, and always an
increase of efficiency when using RQMC instead of MC. However, looking at the
absolute error of the estimate on the right of the Figure, the four different
RQMC estimates show very similar results. This behavior and the fact that the
effective sample size is not a good indicator of the quality of the estimate
in case of quasi-Monte Carlo is not surprising. Computing the ESS requires the
proposals to be independent from one another. This is not the case for RQMC,
as the points are correlated (or even anticorreleted). The comparison between
the iIMC and IS is in favour of IS: while the performances is similar with MC,
IS is much more efficient with RQMC.

For the iMRTH, the results are quite different. The effective sample size
increases with RQMC but the absolute error of the estimate is worse than with
MC. In particular, the mean of the iMRTH chained is heavily biased for the
four RQMC sequences, which is not the case for iIMC. This can be explained as
the two methods present two very different behavior. For an iIMC chain, the
number of replications of a point depends only on the point itself, not on the
past or the future. For a iMRTH chain, the number of replications will depend
on the following proposals from the flow. Therefore, as the RQMC points are
not independent, the whole process is biased. This result is coherent with
previous results on the use of RQMC for Markov chain Monte Carlo [39, 26]. In
these papers, the authors use a different kind of RQMC sequences, called
completely uniformly distributed (CUD) sequences to make RQMC works with MCMC.

##  6 Conclusion

In conclusion, we proposed to replace the classic i.i.d. sample from Monte
Carlo for generating a trained normalizing flow with a randomized quasi-Monte
Carlo sequence. We showed in several experiments that this can lead to a
significant variance reduction and a better convergence rate. The main
drawback is the scaling with the dimension, as the efficiency of RQMC is
rapidly decreasing with the dimension.

The choice of the RQMC sequence is important, and the best choice in our
experiments was a scrambled Sobol’ sequence with a length $2^{k}$. The Box-
Muller transform seems to be better than the inverse transform but this needs
to be confirmed in further experiments. As expected from previous works,
applying RQMC to MCMC methods can lead to biased results, in particular for
the independent Metropolis-Rosenbluth-Teller-Hastings method. The independent
importance markov chain method seems more robust to the use of RQMC, but still
less efficient than importance sampling.

Future works could explore the use RQMC for training normalizing flows and not
only the exploitation of the trained flow. This could be made easier if deep
learning frameworks implement RQMC sequences.

## Acknowledgements

The author was supported by a grant from Région Ile-de-France.

## References

  * [1] M.. Albergo, G. Kanwar and P.. Shanahan  “Flow-Based Generative Models for Markov Chain Monte Carlo in Lattice Field Theory”  In _Physical Review D_ 100.3, 2019, pp. 034515  DOI: 10.1103/PhysRevD.100.034515
  * [2] Charly Andral, Randal Douc, Hugo Marival and Christian P. Robert  “The Importance Markov Chain”  arXiv, 2023  DOI: 10.48550/arXiv.2207.08271
  * [3] George EP Box and Mervin E Muller  “A Note on the Generation of Random Normal Deviates”  In _The annals of mathematical statistics_ 29.2 Institute of Mathematical Statistics, 1958, pp. 610–611
  * [4] R. Cranley and T… Patterson  “Randomization of Number Theoretic Methods for Multiple Integration”  In _SIAM Journal on Numerical Analysis_ 13.6 Society for Industrial and Applied Mathematics, 1976, pp. 904–914  DOI: 10.1137/0713071
  * [5] Josef Dick and Daniel Rudolf  “Discrepancy Estimates for Variance Bounding Markov Chain Quasi-Monte Carlo”  In _Electronic Journal of Probability_ 19, 2014  DOI: 10.1214/EJP.v19-3132
  * [6] Josef Dick, Daniel Rudolf and Houying Zhu  “A Weighted Discrepancy Bound of Quasi-Monte Carlo Importance Sampling”  In _Statistics & Probability Letters_ 149, 2019, pp. 100–106  DOI: 10.1016/j.spl.2019.01.014
  * [7] Josef Dick, Daniel Rudolf and Houying Zhu  “Discrepancy Bounds for Uniformly Ergodic Markov Chain Quasi-Monte Carlo”  In _The Annals of Applied Probability_ 26.5 Institute of Mathematical Statistics, 2016, pp. 3178–3205  DOI: 10.1214/16-AAP1173
  * [8] Laurent Dinh, David Krueger and Yoshua Bengio  “NICE: Non-linear Independent Components Estimation”  arXiv, 2015  arXiv:1410.8516 [cs]
  * [9] Laurent Dinh, Jascha Sohl-Dickstein and Samy Bengio  “Density Estimation Using Real NVP”  arXiv, 2017  arXiv:1605.08803 [cs, stat]
  * [10] Conor Durkan, Artur Bekasov, Iain Murray and George Papamakarios  “Neural Spline Flows”  In _Advances in neural information processing systems_ 32, 2019
  * [11] Marylou Gabrié, Grant M. Rotskoff and Eric Vanden-Eijnden  “Adaptive Monte Carlo Augmented with Normalizing Flows”  In _Proceedings of the National Academy of Sciences_ 119.10, 2022, pp. e2109420119  DOI: 10.1073/pnas.2109420119
  * [12] G Goertzel  “Quota Sampling and Importance Functions in Stochastic Solution of Particle Problems”, 1949  URL: https://www.osti.gov/biblio/4405957
  * [13] John A. Gregory and Roger Delbourgo  “Piecewise Rational Quadratic Interpolation to Monotonic Data”  In _IMA Journal of Numerical Analysis_ 2.2 Oxford University Press, 1982, pp. 123–130
  * [14] Louis Grenioux, Alain Durmus, Éric Moulines and Marylou Gabrié  “On Sampling with Approximate Transport Maps”  arXiv, 2023  arXiv:2302.04763 [cs, stat]
  * [15] John H Halton  “On the Efficiency of Certain Quasi-Random Sequences of Points in Evaluating Multi-Dimensional Integrals”  In _Numerische Mathematik_ 2 Springer, 1960, pp. 84–90
  * [16] W.. Hastings  “Monte Carlo Sampling Methods Using Markov Chains and Their Applications”  In _Biometrika_ 57.1, 1970, pp. 97–109  DOI: 10.1093/biomet/57.1.97
  * [17] Zhijian He, Zhan Zheng and Xiaoqun Wang  “On the Error Rate of Importance Sampling with Randomized Quasi-Monte Carlo”  In _SIAM Journal on Numerical Analysis_ 61.2, 2023, pp. 515–538  DOI: 10.1137/22M1510121
  * [18] Edmund Hlawka  “Funktionen von beschränkter variatiou in der theorie der gleichverteilung”  In _Annali di Matematica Pura ed Applicata_ 54.1 Springer, 1961, pp. 325–333
  * [19] Li Huang and Lei Wang  “Accelerated Monte Carlo Simulations with Restricted Boltzmann Machines”  In _Physical Review B_ 95.3 American Physical Society, 2017, pp. 035105  DOI: 10.1103/PhysRevB.95.035105
  * [20] Stephen Joe and Frances Y. Kuo  “Constructing Sobol Sequences with Better Two-Dimensional Projections”  In _SIAM Journal on Scientific Computing_ 30.5, 2008, pp. 2635–2654  DOI: 10.1137/070709359
  * [21] Herman Kahn  “Stochastic (Monte Carlo) Attenuation Analysis. June 14, 1949. R- 163.”, 1949
  * [22] Xiongwen Ke et al.  “Adaptive Importance Sampling and Quasi-Monte Carlo Methods for 6G URLLC Systems”  arXiv, 2023  arXiv:2303.03575 [eess, stat]
  * [23] Ivan Kobyzev, Simon J.D. Prince and Marcus A. Brubaker  “Normalizing Flows: An Introduction and Review of Current Methods”  In _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 43.11, 2021, pp. 3964–3979  DOI: 10.1109/TPAMI.2020.2992934
  * [24] JF Koksma  “Een algemeene stelling uit de theorie der gelijkmatige verdeeling modulo 1”  In _Mathematica B (Zutphen)_ 11.7-11, 1942, pp. 43
  * [25] Pierre l’Ecuyer, Marvin K Nakayama, Art B Owen and Bruno Tuffin  “Confidence Intervals for Randomized Quasi-Monte Carlo Estimators”  In _Winter Simulation Conference_ , 2023
  * [26] Sifan Liu  “Langevin Quasi-Monte Carlo”  arXiv, 2023  DOI: 10.48550/arXiv.2309.12664
  * [27] Jiří Matoušek  “On the L2-Discrepancy for Anchored Boxes”  In _Journal of Complexity_ 14.4, 1998, pp. 527–556  DOI: 10.1006/jcom.1998.0489
  * [28] Nicholas Metropolis et al.  “Equation of State Calculations by Fast Computing Machines”  In _The Journal of Chemical Physics_ 21.6 American Institute of Physics, 1953, pp. 1087–1092  DOI: 10.1063/1.1699114
  * [29] Nicholas Metropolis and S. Ulam  “The Monte Carlo Method”  In _Journal of the American Statistical Association_ 44.247 [American Statistical Association, Taylor & Francis, Ltd.], 1949, pp. 335–341  JSTOR: 2280232
  * [30] Laurence Illing Midgley et al.  “Flow Annealed Importance Sampling Bootstrap”  In _The Eleventh International Conference on Learning Representations_ , 2023
  * [31] Thomas Müller et al.  “Neural Importance Sampling”  In _ACM Transactions on Graphics_ 38.5, 2019, pp. 145:1–145:19  DOI: 10.1145/3341156
  * [32] Frank Noé, Simon Olsson, Jonas Köhler and Hao Wu  “Boltzmann Generators: Sampling Equilibrium States of Many-Body Systems with Deep Learning”  In _Science_ 365.6457, 2019, pp. eaaw1147  DOI: 10.1126/science.aaw1147
  * [33] Art B. Owen  “A Randomized Halton Algorithm in R”  arXiv, 2017  DOI: 10.48550/arXiv.1706.02808
  * [34] Art B. Owen  “Halton Sequences Avoid the Origin”  In _SIAM Review_ 48.3, 2006, pp. 487–503  DOI: 10.1137/S0036144504441573
  * [35] Art B. Owen  “Monte Carlo Theory, Methods and Examples”, 2013
  * [36] Art B. Owen  “On Dropping the First Sobol’ Point”  In _Monte Carlo and Quasi-Monte Carlo Methods_ , Springer Proceedings in Mathematics & Statistics  Cham: Springer International Publishing, 2022, pp. 71–86  DOI: 10.1007/978-3-030-98319-2˙4
  * [37] Art B. Owen  “Randomly Permuted $(t,m,s)$-Nets and $(t,s)$-Sequences”  In _Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing_ New York: Springer-Verlag, 1995, pp. 299–317
  * [38] Art B. Owen  “Scrambled Net Variance for Integrals of Smooth Functions”  In _The Annals of Statistics_ 25.4 Institute of Mathematical Statistics, 1997, pp. 1541–1562
  * [39] Art B. Owen and Seth D. Tribble  “A Quasi-Monte Carlo Metropolis Algorithm”  In _Proceedings of the National Academy of Sciences_ 102.25, 2005, pp. 8844–8849  DOI: 10.1073/pnas.0409596102
  * [40] George Papamakarios et al.  “Normalizing Flows for Probabilistic Modeling and Inference”  In _The Journal of Machine Learning Research_ 22.1, 2021, pp. 57:2617–57:2680
  * [41] Matthew D. Parno and Youssef M. Marzouk  “Transport Map Accelerated Markov Chain Monte Carlo”  In _SIAM/ASA Journal on Uncertainty Quantification_ 6.2, 2018, pp. 645–682  DOI: 10.1137/17M1134640
  * [42] Christian P. Robert and George Casella  “Monte Carlo Statistical Methods”, Springer Texts in Statistics  New York, NY: Springer New York, 2010  DOI: 10.1007/978-1-4757-4145-2
  * [43] Pamphile T. Roy, Art B. Owen, Maximilian Balandat and Matt Haberland  “Quasi-Monte Carlo Methods in Python”  In _Journal of Open Source Software_ 8.84, 2023, pp. 5309  DOI: 10.21105/joss.05309
  * [44] Il’ya Meerovich Sobol’  “On the Distribution of Points in a Cube and the Approximate Evaluation of Integrals”  In _Zhurnal Vychislitel’noi Matematiki i Matematicheskoi Fiziki_ 7.4 Russian Academy of Sciences, Branch of Mathematical Sciences, 1967, pp. 784–802
  * [45] E.. Tabak and Cristina V. Turner  “A Family of Nonparametric Density Estimation Algorithms”  In _Communications on Pure and Applied Mathematics_ 66.2, 2013, pp. 145–164  DOI: 10.1002/cpa.21423
  * [46] Esteban G. Tabak and Eric Vanden-Eijnden  “Density Estimation by Dual Ascent of the Log-Likelihood”  In _Communications in Mathematical Sciences_ 8.1 International Press of Boston, 2010, pp. 217–233
  * [47] Florian Wenzel, Alexander Buchholz and Stephan Mandt  “Quasi-Monte Carlo Flows”  In _Proceedings of the 3rd Workshop on Bayesian Deep Learning_ , 2018
  * [48] Kaze WK Wong, Marylou Gabrié and Daniel Foreman-Mackey  “flowMC: Normalizing-flow Enhanced Sampling Package for Probabilistic Inference in Jax”  In _arXiv preprint arXiv:2211.06397_ , 2022  arXiv:2211.06397