# Retrieval Based Time Series Forecasting

Baoyu Jing<EMAIL_ADDRESS>University of Illinois at Urbana-Champaign ,
Si Zhang<EMAIL_ADDRESS>Meta ,  Yada Zhu<EMAIL_ADDRESS>IBM Research ,  Bin
Peng<EMAIL_ADDRESS>University of Illinois at Urbana-Champaign ,  Kaiyu
Guan<EMAIL_ADDRESS>University of Illinois at Urbana-Champaign ,  Andrew
Margenot<EMAIL_ADDRESS>University of Illinois at Urbana-Champaign and
Hanghang Tong<EMAIL_ADDRESS>University of Illinois at Urbana-Champaign

(2022)

###### Abstract.

Time series data appears in a variety of applications such as smart
transportation and environmental monitoring. One of the fundamental problems
for time series analysis is time series forecasting. Despite the success of
recent deep time series forecasting methods, they require sufficient
observation of historical values to make accurate forecasting. In other words,
the ratio of the output length (or forecasting horizon) to the sum of the
input and output lengths should be low enough (e.g., 0.3). As the ratio
increases (e.g., to 0.8), the uncertainty for the forecasting accuracy
increases significantly. In this paper, we show both theoretically and
empirically that the uncertainty could be effectively reduced by retrieving
relevant time series as references. In the theoretical analysis, we first
quantify the uncertainty and show its connections to the Mean Squared Error
(MSE). Then we prove that models with references are easier to learn than
models without references since the retrieved references could reduce the
uncertainty. To empirically demonstrate the effectiveness of the retrieval
based time series forecasting models, we introduce a simple yet effective two-
stage method, called ReTime consisting of a relational retrieval and a content
synthesis. We also show that ReTime can be easily adapted to the spatial-
temporal time series and time series imputation settings. Finally, we evaluate
ReTime on real-world datasets to demonstrate its effectiveness.

††copyright: acmcopyright††journalyear: 2022††doi:
XXXXXXX.XXXXXXX††conference: CIKM workshop on Applied Machine Learning Methods
for Time Series Forecasting; October 17–21, 2022; Atlanta, GA,
USA.††booktitle: CIKM’22 workshop on Applied Machine Learning Methods for Time
Series Forecasting (AMLTS), October 17–21, 2022, Atlanta, GA, USA††price:
15.00††isbn: 978-1-4503-XXXX-X/18/06

##  1\. Introduction

Time series analysis has received paramount interest in numerous real-world
applications (Zhou et al., 2021; Gamboa, 2017; Li et al., 2019; Wu et al.,
2021; Zhou et al., 2020; Zhu et al., 2022; Yang et al., 2019), such as smart
transportation and environmental monitoring. Accurately forecasting time
series provides valuable insights for making transport policies (Li et al.,
2018a) and climate change policies (Sarker et al., 2012).

One of the fundamental problems for time series analysis is time series
forecasting. Despite the success of recent deep learning methods for time
series forecasting (Zhou et al., 2021; Gamboa, 2017; Li et al., 2019; Wu et
al., 2021; Zhou et al., 2020), a sufficient observation of the input
historical time series is required. Specifically, the ratio of the output
length (i.e., forecasting horizon) to the sum of the input and output length
should be sufficiently low (e.g., 0.3). In fact, this ratio is a special case
of the missing rate used in time series imputation, and thus for clarity and
consistency with the literature, we formulate the task of time series
forecasting from the perspective of time series imputation in this paper.
Please refer to Section 2 for details. When the missing rate increases to a
high level (e.g., 0.8), the uncertainty of the forecasting accuracy will
increase significantly. In real-world applications, it is common that one
wants to forecast future values based on very limited observations. In smart
transportation, government administrators might be interested in forecasting
the traffic conditions of a road with broken sensors (Wang and Kockelman,
2009). In environmental monitoring, geologists always have a desire of
obtaining the temperature of a specific location, where sensors cannot be
easily placed (Nalder and Wein, 1998).

To address this problem, we first formally quantify the uncertainty of the
forecasting based on the entropy of the ground truth conditioned on the
predicted values, and show its connections to the Mean Squared Error (MSE).
Then we theoretically prove that models with reference time series are easier
to train than those without references, since the references could reduce
uncertainties. Motivated by the theoretical analysis, we introduce a simple
yet effective two-stage time forecasting method called ReTime. Given a target
time series, in the first relational retrieval stage, ReTime retrieves the
references from a database based on the relations among time series. We use
relational retrieval rather than content based retrieval since the input
historical values of the target time series could be very unreliable when the
input length is very short. Thus, content-based methods could retrieve
unreliable references. In comparison, the relational information is usually
reliable and easy to obtain in practice (Bonnet et al., 2001; Pelkonen et al.,
2015; Rhea et al., 2017; Song et al., 2018a), such as whether two sensors are
adjacent to each other in traffic/environmental monitoring. In the second
content synthesis stage, ReTime synthesizes the future values based on the
content of the target and the references. Next, we show that the proposed
ReTime could be easily applied to the time series imputation task and the
spatial-temporal time series setting. Finally, we empirically evaluate ReTime
on two real-world datasets to demonstrate its effectiveness.

The main contributions of the paper are summarized as follows:

  * •

Theoretical Analysis. We theoretically quantify the uncertainty of the
predicted values based on conditional entropy and show its connections to the
MSE. We also theoretically demonstrate that models with references are easier
to train than those without references.

  * •

Algorithm. We introduce a two-stage method ReTime for time series forecasting,
which is comprised of relational retrieval and content synthesis. ReTime can
also be easily applied to the spatial-temporal time series and time series
imputation settings.

  * •

Empirical Evaluation. We evaluate ReTime on two real-world datasets and
various settings (i.e., single/spatial temporal forecasting/imputation) to
demonstrate its effectiveness.

##  2\. Preliminary

The ratio ${L_{out}}/{(L_{in}+L_{out})}$ of the output length $L_{out}$ to the
sum of the input length $L_{in}$ and output length $L_{out}$ is a special case
of the missing rate used in time series imputation, and thus we formulate
tasks of time series forecasting from the perspective of time series
imputation. We summarize the mathematical notations in Table 1.

###### Definition 2.1 (Single Time Series Forecasting).

Given an incomplete target time series $\mathbf{X}\in\mathbb{R}^{T\times v}$,
where $T$ and $v$ are the numbers of time steps and variates, along with its
indicator mask $\mathbf{M}\in\\{0,1\\}^{{T}\times v}$, which indicates the
absence/presence of the data points, the task aims to generate a new target
time series $\hat{\mathbf{X}}\in\mathbb{R}^{T\times v}$ to predict the missing
values $(1-\mathbf{M})\odot\hat{\mathbf{X}}$.

###### Definition 2.2 (Retrieval Based Single Time Series Forecasting).

Given an incomplete target time series $\mathbf{X}\in\mathbb{R}^{T\times v}$,
where $T$ and $v$ are the numbers of time steps and variates, along with its
indicator mask $\mathbf{M}\in\\{0,1\\}^{{T}\times v}$, the task aims to
generate a new target time series $\hat{\mathbf{X}}\in\mathbb{R}^{T\times v}$
to predict the missing values $(1-\mathbf{M})\odot\hat{\mathbf{X}}$ based on
the input target time series $\mathbf{X}$ and the $K$ reference time series
$\\{\mathbf{Y}_{k}\\}_{k=1}^{K}$ retrieved from a database
$\\{\mathbf{Y}_{n}^{\prime}\in\mathbb{R}^{T^{\prime}\times v}\\}_{n=1}^{N}$,
where $N\gg K$ is the number of time series and $T^{\prime}\gg T$.

###### Definition 2.3 (Spatial-Temporal Time Series Forecasting).

Given an incomplete spatial-temporal time series
$\mathbf{X}\in\mathbb{R}^{N\times T\times v}$, where $N$, $T$ and $v$ are the
numbers of time series, time steps, and variates, along with its indicator
mask $\mathbf{M}\in\\{0,1\\}^{N\times T\times v}$, where 0/1 indicates the
absence/presence of the data points, and its adjacency matrix
$\mathbf{A}\in\mathbb{R}^{N\times N}$, the task is to generate a new time
series $\hat{\mathbf{X}}\in\mathbb{R}^{N\times T\times v}$ to predict the
missing values $(1-\mathbf{M})\odot\hat{\mathbf{X}}$.

Note that for forecasting, missing points and zeros are concentrated at the
end of the target $\mathbf{X}$ and mask $\mathbf{M}$ after a certain
separation time step $\tau$: $\mathbf{M}[\tau:T]=0$.

Notation | Description
---|---
$\mathbf{X}$ | incomplete target time series
$\tilde{\mathbf{X}}$ | complete target time series
$\hat{\mathbf{X}}$ | predicted target time series
$\mathbf{Y}_{k}$ |  $k$-th reference time series
$\mathbf{Y}_{n}^{\prime}$ |  $n$-th time series in the database
$\mathbf{M}$ | indicator mask of the target time series
$\mathbf{A}^{\prime}$ | adjacency matrix for the database $\\{\mathbf{Y}_{n}^{\prime}\\}_{n=1}^{N}$
$\mathbf{A}$ | adjacency matrix for $\\{\mathbf{Y}_{n}^{\prime}\\}_{n=1}^{N}$ and $\mathbf{X}$
$\mathcal{R}$ | relation between $\mathbf{X}$ and $\\{\mathbf{Y}_{n}^{\prime}\\}_{n=1}^{N}$
$T$ | length of time series
$T^{\prime}$ | length of time series in $\\{\mathbf{Y}_{n}^{\prime}\\}_{n=1}^{N}$
$N$ | number of time series in the database
$K$ | number of reference time series
$v$ | number of variates
$d$ | size of hidden dimension
$\tau$ | separation time separating history and future time
Table 1. Mathematical Notations

##  3\. Theoretical Motivation

(a) With references

(b) Without references

Figure 1. Graphical models for methods w. or w/o references. $X$, $\hat{X}$,
$\tilde{X}$ and $Y$ are random variables for incomplete, generated, complete,
and retrieved time series. $g$ and $r$ are the generation and retrieval
models. $p$ is the relation between $\hat{X}$ and $\tilde{X}$. Figure 2.
Overview of ReTime. Given $\mathbf{X}$ and
$\\{\mathbf{Y}_{n}^{\prime}\\}_{n=1}^{N}$, ReTime first retrieves the top $K$
references $\\{\mathbf{Y}_{k}\\}_{k=1}^{K}$ based on the relations
$\mathbf{A}^{\prime}$ and $\mathcal{R}$, and then combines the content of
$\mathbf{X}$ and $\\{\mathbf{Y}_{k}\\}_{k=1}^{K}$ to generate
$\hat{\mathbf{X}}$. Solid/dashed curves are observed/unobserved values.

An illustration of graphical models for methods with or without references is
presented in Figure 1, which shows the relations among random variables. $X$,
$\tilde{X}$, $\hat{X}$ and $Y$ denote the random variables for the incomplete
target, complete target, generated target, and retrieved reference time series
respectively. $g$ and $r$ denote the generation and retrieval model
respectively. $p$ is the relation between $\hat{X}$ and $\tilde{X}$. Following
the common practice for linear regression (Murphy, 2012), which assumes a
Gaussian noise between $\hat{X}$ and $\tilde{X}$, we define the relation $p$
as:

(1)
$p(\tilde{\mathbf{x}}|\hat{\mathbf{x}})=\mathcal{N}(\tilde{\mathbf{x}}|\hat{\mathbf{x}},\sigma^{2}\mathbf{I})$

where $\mathcal{N}$ denotes the Gaussian distribution, $\hat{\mathbf{x}}$ and
$\tilde{\mathbf{x}}\in\mathbb{R}^{v}$ denote the values of $\hat{\mathbf{X}}$
and $\tilde{\mathbf{X}}\in\mathbb{R}^{(T-\tau)\times v}$ at a future time step
$t\in[\tau,T]$, $\sigma$ is the standard deviation, and
$\mathbf{I}\in\mathbb{R}^{v\times v}$ is the identity matrix. Note that
$\mathbf{X}[:\tau]=\tilde{\mathbf{X}}[:\tau]$ for historical steps
$t\in[1,\tau)$, and thus we only consider the future time steps
$t\in[\tau,T]$.

We first quantify the uncertainty $\Delta$ for the accuracy of $\hat{X}$ in
Definition 3.1 as the entropy of $\tilde{X}$ conditioned on $\hat{X}$.
According to Equation 1 and the definition of conditional entropy, we can
calculate the uncertainty $\Delta$ (Lemma 3.1). As the only parameter in
$\Delta$ is the standard deviation $\sigma$ in Equation 1, we can further
prove that $\Delta$ is equivalent to the MSE between $\tilde{X}$ and
$\hat{X}$.

###### Definition 3.1 (Uncertainty of $\hat{X}$).

(2) $\Delta=H(\tilde{X}|\hat{X})$

where $H$ denotes the entropy.

###### Lemma 3.1 (Uncertainty Calculation).

According to the definition of conditional entropy and Equation (1), we have:

(3) $\Delta=\frac{v}{2}(1+\log 2\pi\sigma^{2})$

###### Lemma 3.2 (Equivalence between Uncertainty and MSE).

The uncertainty is equivalent to MSE.

(4) $\Delta\Leftrightarrow MSE$

###### Proof.

The only parameter of $\Delta$ in Lemma 3.1 is the standard deviation
$\sigma$, which can be estimated by:

(5)
$\sigma=\sqrt{\frac{1}{Z}\sum_{z=1}^{Z}||\tilde{\mathbf{x}}_{z}-\hat{\mathbf{x}}_{z}||^{2}}$

where $Z$ is the total number of data pairs. The item under the square root is
MSE between $\tilde{X}$ and $\hat{X}$. ∎

We further study the relations among the inputs $X$, $Y$, output $\hat{X}$ of
the generation model $g$, and the complete time series $\tilde{X}$ based on
their dependencies shown in Figure 1. Firstly, given $\tilde{X}$ and
$\hat{X}$, we show in Lemma 3.3 that minimizing their MSE loss is equivalent
to maximizing their mutual information $I(\tilde{X};\hat{X})$. Secondly, we
prove that adding the retrieved reference $Y$ to the input could reduce the
uncertainty for $X$. Lemma 3.4 shows that methods with $Y$ have a higher
lower-bound than those without $Y$ for the mutual information of the ground-
truth $\tilde{X}$ and the predicted values $\hat{X}$:
$I(\tilde{X};\hat{X})\geq I(\tilde{X};X,Y)\geq I(\tilde{X};X)$. Finally, due
to the equivalence of MSE and MI (Lemma 3.3), we can conclude that models with
$Y$ (Figure 1(a)) are easier to learn than models without $Y$ (Figure 1(b))
under the MSE loss.

###### Lemma 3.3 (Equivalence between MSE and MI).

Minimizing the MSE loss of $\tilde{X}$ and $\hat{X}$ is equivalent to maximize
the mutual information of $\tilde{X}$ and $\hat{X}$: $I(\tilde{X};\hat{X})$.

(6)
$\min\mathbb{E}_{p(\tilde{\mathbf{x}},\hat{\mathbf{x}})}||\tilde{\mathbf{x}}-\hat{\mathbf{x}}||^{2}\Leftrightarrow\max
I(\tilde{X},\hat{X})$

where $p(\tilde{\mathbf{x}},\hat{\mathbf{x}})$ indicates whether
$\tilde{\mathbf{x}}$ and $\hat{\mathbf{x}}$ is a true pair.

###### Proof.

Minimizing MSE of $\tilde{\mathbf{x}}$ and $\hat{\mathbf{x}}$ is equivalent to
maximizing the log-likelihood $\log p(\tilde{\mathbf{x}}|\hat{\mathbf{x}})$
(Murphy, 2012), where $p$ is given in Equation (1):

(7)
$\min\mathbb{E}_{p(\tilde{\mathbf{x}},\hat{\mathbf{x}})}||\tilde{\mathbf{x}}-\hat{\mathbf{x}}||^{2}\Leftrightarrow\max\mathbb{E}_{p(\tilde{\mathbf{x}},\hat{\mathbf{x}})}[\log
p(\tilde{\mathbf{x}}|\hat{\mathbf{x}})]$

Therefore, we only need to prove

(8) $\max\mathbb{E}_{p(\tilde{\mathbf{x}},\hat{\mathbf{x}})}[\log
p(\tilde{\mathbf{x}}|\hat{\mathbf{x}})]\Leftrightarrow\max
I(\tilde{X};\hat{X})$

In fact,

(9) $I(\tilde{X};\hat{X})=H(\tilde{X})-H(\tilde{X}|\hat{X})$

and $H(\tilde{X})$ is a constant since the ground-truth $\tilde{X}$ is fixed
in the dataset. Thus, we have

(10) $\max I(\tilde{X};\hat{X})\Leftrightarrow\max-H(\tilde{X}|\hat{X})$

According to the definition of conditional entropy, we have

(11)
$-H(\tilde{X}|\hat{X})=\mathbb{E}_{p(\tilde{\mathbf{x}},\hat{\mathbf{x}})}[\log
p(\tilde{\mathbf{x}}|\hat{\mathbf{x}})]$

The proof is concluded by combining Equations (7)(10)(11). ∎

###### Lemma 3.4 (MI Monotonicity).

The following MI inequalities hold for the graphical model shown in Figure
1(a).

(12) $I(\tilde{X};\hat{X})\geq I(\tilde{X};X,Y)\geq I(\tilde{X};X)$

###### Sketch of Proof.

The first inequality is derived based on the data processing inequality
(Mezard and Montanari, 2009). The second inequality holds since
$I(\tilde{X};X,Y)=I(\tilde{X};X)+I(\tilde{X};Y|X)$ and $I(\tilde{X};Y|X)\geq
0$. ∎

##  4\. Methodology

Figure 2 is an overview of ReTime. In the first stage, ReTime retrieves the
top $K$ references $\\{\mathbf{Y}_{k}\\}_{k=1}^{K}$ for the target
$\mathbf{X}$ from the database $\\{\mathbf{Y}_{n}^{\prime}\\}_{n=1}^{N}$,
based on the relations $\mathcal{R}$ between the target and database and the
relation graph $\mathbf{A}^{\prime}$ of the database. In the second stage,
ReTime combines $\mathbf{X}$ and $\\{\mathbf{Y}_{k}\\}_{k=1}^{K}$ to generate
$\hat{\mathbf{X}}$.

###  4.1. Relational Retrieval

When the missing rate of the target $\mathbf{X}$ is high (e.g., 0.8), it will
be hard to accurately complete $\mathbf{X}$ merely based on the observed
historical content of $\mathbf{X}$. Even for recent forecasting methods (Zhou
et al., 2021), the uncertainty of the prediction accuracy could be very high
under a high missing rate. To reduce the uncertainty, we propose to retrieve
$K$ references $\\{\mathbf{Y}_{k}\\}_{k=1}^{K}$ from the database
$\\{\mathbf{Y}_{n}^{\prime}\\}_{n=1}^{N}$, based on the relations
$\mathcal{R}$ and $\mathbf{A}^{\prime}$. We choose relational retrieval over
content-based retrieval since the observed historical content could be noisy
and unreliable when the missing rate is high.

Given $\mathbf{A}^{\prime}\in\mathbb{R}^{N\times N}$ and $\mathcal{R}$, we
first construct a new adjacency matrix
$\mathbf{A}\in\mathbb{R}^{(N+1)\times(N+1)}$ by appending $\mathcal{R}$ to the
last row and column of $\mathbf{A}^{\prime}$. Then we use the Random Walk with
Restart (RWR) (Tong et al., 2006) to obtain the relational proximity scores
between $\mathbf{X}$ and
$\\{\mathbf{Y}_{n}^{\prime}\in\mathbb{R}^{T^{\prime}\times v}\\}_{n=1}^{N}$,
whose closed-form solution is given by:

(13) $\mathbf{p}=(1-c)(\mathbf{I}-c\tilde{\mathbf{A}})^{-1}\mathbf{e}$

where $\tilde{\mathbf{A}}$ is the normalized adjacency matrix, $\mathbf{I}$ is
the identity matrix and $c\in(0,1)$ is the tunable damping factor;
$\mathbf{e}\in\\{0,1\\}^{N+1}$ is the indicator vector of $\mathbf{X}$, where
$\mathbf{e}[N+1]=1$ and $\mathbf{e}[n]=0$ for $\forall n\in[1,\cdots,N]$;
$\mathbf{p}\in\mathbb{R}^{N+1}$ is the relational proximity vector where
$\mathbf{p}[n]$ describes the proximity between $\mathbf{X}$ and
$\mathbf{Y}_{n}^{\prime}$. Given $\mathbf{p}$, we retrieve the top $K$ time
series $\\{\mathbf{Y}_{k}\\}_{k=1}^{K}$ as the references.

Instead of using the entire $\mathbf{Y}_{k}\in\mathbb{R}^{T^{\prime}\times v}$
as the reference, we only use a $T$-length ($T\ll T^{\prime}$) snippet of it,
since using the entire sequence could introduce much irrelevant noisy
information. Let $T_{S}$ and $T_{S}^{\prime}$ be the start time of the target
$\mathbf{X}$ and the reference snippet respectively. Many time series have
clear periodicity patterns, such as the weekly pattern of traffic monitoring
data and the yearly pattern of the air temperature. In this paper, we exploit
the periodicity patterns and set the time difference between the target and
$\Delta T=T_{S}-T_{S}^{\prime}$ as the length of one period.

###  4.2. Content Synthesis

Despite the relational closeness of $\mathbf{X}$ and
$\\{\mathbf{Y}_{k}\\}_{k=1}^{K}$, they usually have content discrepancies. For
example, peaks and valleys of $\mathbf{X}$ and $\mathbf{Y}_{k}$ might be
different, and $\mathbf{Y}_{k}$ might be noisy. Besides, there are also
temporal dependencies among different time steps. Therefore, it is necessary
to build a model to combine their content.

Figure 3 presents an illustration of the content synthesis model, which is
comprised of the input, aggregation, and output modules. The input module maps
$\mathbf{X}\in\mathbb{R}^{T\times v}$ and
$\\{\mathbf{Y}_{k}\in\mathbb{R}^{T\times v}\\}_{k=1}^{K}$ into embeddings
$\mathbf{H}\in\mathbb{R}^{(K+1)\times T\times d}$, where $d$ is the size of
hidden dimension. Then the aggregation module aggregates the content of
$\mathbf{H}$ across the $K+1$ time series and $T$ time steps into the
aggregated embeddings $\mathbf{H}$. Finally, the output module generates the
completed time series $\hat{\mathbf{X}}$ based on the aggregated embeddings
$\mathbf{H}$.

Figure 3. Content Synthesis

A - Input Module. The input module maps the time series into the embedding
space. Given $\mathbf{X}$ and $\\{\mathbf{Y}_{k}\\}_{k=1}^{K}$, we first apply
separate linear layers to them. Then we apply the position encoding (Vaswani
et al., 2017) and layer norm (Ba et al., 2016) to obtain the embeddings
$\mathbf{H}\in\mathbb{R}^{T\times d}$ and
$\\{\mathbf{H}^{\prime}_{k}\in\mathbb{R}^{T\times d}\\}_{k=1}^{K}$. Finally,
they are concatenated into $\hat{\mathbf{H}}\in\mathbb{R}^{(K+1)\times T\times
d}$, where
$\hat{\mathbf{H}}[1:K]=[\mathbf{H}^{\prime}_{1},\dots,\mathbf{H}_{K}^{\prime}]$
and $\hat{\mathbf{H}}[K+1]={\mathbf{H}}$.

B - Aggregation Module. The aggregation module jointly considers the content
discrepancies between $\mathbf{H}$ and
$\\{\mathbf{H}^{\prime}_{k}\\}_{k=1}^{K}$ at each time step and the temporal
dependencies across time steps. Based on the multi-head self-attention
(Vaswani et al., 2017), we build content and temporal attention models to
handle the content discrepancies and temporal dependencies. The aggregation
module sequentially applies the content and temporal attention models to
$\hat{\mathbf{H}}\in\mathbb{R}^{(K+1)\times T\times d}$, which is comprised of
$L$ blocks, and the structure of the block is shown in Figure 3. Within the
block, firstly, the content attention computes content attention scores over
the first dimension ($K+1$ time series) of $\hat{\mathbf{H}}$ for each step
$t\in[1,\cdots,T]$, and produces the content embeddings
$\mathbf{Z}\in\mathbb{R}^{(K+1)\times T\times d}$. Secondly, the temporal
attention computes temporal attention scores over the second dimension ($T$
steps) of $\mathbf{Z}$ for each time series $k\in[1,\cdots,K+1]$, and encodes
temporal information into $\mathbf{Z}$. Finally, a forward layer (Vaswani et
al., 2017) is applied to $\mathbf{Z}$ to obtain the aggregated embeddings
$\hat{\mathbf{H}}\in\mathbb{R}^{(K+1)\times T\times d}$.

C - Output Module. The output module maps the aggregated embeddings
$\hat{\mathbf{H}}$ from the embedding space $\mathbb{R}^{d}$ back into the
original space $\mathbb{R}^{V}$ of the time series. Specifically, the output
module takes input as $\hat{\mathbf{H}}[K+1]$, which is the embeddings of the
target $\mathbf{X}$. Then a Multi-Layer Perceptron (MLP) is applied over
$\mathbf{H}[K+1]$ to generate the predicted target time series
$\hat{\mathbf{X}}$.

###  4.3. Adaptation to Other Settings

A - Spatial-Temporal Time Series. Spatial-temporal time series is ubiquitous
and have attracted a lot of attention (Yu et al., 2016; Cao et al., 2020; Wu
et al., 2020; Jing et al., 2021a). Different from the retrieval-based single
time series forecasting, where the target is a single time series _out of_ the
database $\mathbf{X}\notin\\{\mathbf{Y}_{n}\\}_{n=1}^{N}$, the target of the
spatial-temporal time series is the entire dataset
$\mathbf{X}=\\{\mathbf{Y}_{n}\\}_{n=1}^{N}$. ReTime can be naturally adapted
to spatial-temporal time series by treating each single time series of
spatial-temporal time series as the target and the rest time series as the
database. The relational retrieval can be interpreted as the diffusion graph
kernel (Tong et al., 2006), which identifies the most important neighbors, and
the content synthesis aggregates the information of the neighbors.

B - Time Series Imputation. It is obvious that ReTime can be naturally applied
to the time series imputation task for regularly-sampled time series.

###  4.4. Training

During training, given a complete target time series $\tilde{\mathbf{X}}$, we
generate a binary mask $\mathbf{M}$ to obtain the incomplete target time
series $\mathbf{X}=\mathbf{M}\odot\tilde{\mathbf{X}}$, where $\odot$ is the
Hadamard product. For forecasting, the values after the pre-defined time step
$\tau$ are set as zeros. For imputation, we randomly generate the values for
the masks according to the pre-defined missing rate. Then we feed $\mathbf{X}$
with its references $\\{\mathbf{Y}_{k}\\}_{k=1}^{K}$, which are retrieved by
the relational retrieval model, to the content synthesis model to obtain
$\hat{\mathbf{X}}$. We use the standard Mean Squared Error (MSE) between
$\hat{\mathbf{X}}$ and $\tilde{\mathbf{X}}$ to train the model.

(a) Traffic Forecasting

(b) Temperature Forecasting

(c) Traffic Imputation

(d) Temperature Imputation

Figure 4. RMSE scores on single time series forecasting and imputation. The
lower the better.

(a) Traffic Forecasting

(b) Temperature Forecasting

(c) Traffic Imputation

(d) Temperature Imputation

Figure 5. RMSE scores on spatial-temporal time series forecasting and
imputation. The lower the better.

##  5\. Experiments

###  5.1. Experimental Setup

Datasets. We evaluate ReTime on two real-world datasets, where the shapes of
the datasets are formulated as $N\times T\times v$:

Traffic dataset is collected from Caltrans Performance Measurement System
(PeMS)111https://dot.ca.gov/programs/traffic-operations/mpr/pems-source. It
contains hourly average speed and occupancy collected from 2,000 sensor
stations in District 7 of California during June 2018. The size of the dataset
is $2000\times 720\times 2$. The relation between two stations is whether they
are adjacent.

Temperature is a subset of version 3 of the 20th Century
Reanalysis222https://psl.noaa.gov/data/gridded/data.20thC_ReanV3.html data
(Slivinski et al., 2019). It contains the monthly average temperature from
2001 to 2015, which covers a $30\times 30$ area of North America, from
$30^{\circ}$ N to $60^{\circ}$ N, $80^{\circ}$ W to $110^{\circ}$ W. The shape
is $900\times 180\times 1$. The relation between two locations is whether they
are adjacent.

For the single time series setting, given $N$ time series, we randomly select
10%/10% time series for validation and testing. As a result, the
train/validation/test splits are 1600/200/200 and 720/90/90 for the traffic
and temperature datasets. For the spatial-temporal setting, we select 10%/10%
snippets from the entire datasets for validation/test. After splitting data,
we segment each time series into 24/12-length snippets for traffic/temperature
datasets respectively, which cover one day/one year. Time series are
normalized according to the mean and standard deviation of the training sets.

Evaluation Tasks. We evaluate ReTime for both single time series and spatial-
temporal time series, and on both forecasting and imputation tasks. For each
setting, we fix the snippet length of the target and change the missing rate
from 0.2 to 0.8.

Comparison Methods. We compare ReTime with the following baselines. The
methods for single time series setting include block-style methods: N-BEATS
(Oreshkin et al., 2020) and Informer (Zhou et al., 2021), and RNN based
methods: BRITS (Cao et al., 2018a), E2GAN (Luo et al., 2019). The methods for
spatial-temporal time series setting include block-style methods: StemGNN (Cao
et al., 2020) and MTGNN (Wu et al., 2020), and RNN based methods: DCRNN (Li et
al., 2018b), NET3 (Jing et al., 2021a). Additionally, we also use the
reference (REF) with the highest relational score to the target as a baseline.

Implementation Details. The hidden dimensions for the traffic/temperature
datasets are 256/64. The numbers of blocks and attention heads are 8/4 for the
traffic/temperature datasets. The learning rates are tuned within [0.001,
0.0001]. Early stopping is applied on the validation set to prevent over-
fitting. $K$ is tuned within [1, 5, 10, 20]. For the relational retrieval, we
set the damping factor of RWR as $c=0.9$. Given a target time series snippet,
which starts from the time $T_{S}$, we select reference snippets starting from
$T_{S}^{\prime}$, and denote their difference as $\Delta
T=T_{S}-T_{S}^{\prime}$. For forecasting, $\Delta T$ is one week/one year for
the traffic/temperature datasets, which is generally the length of one cycle.
For imputation, we set $\Delta T=0$, which has the best performance in
experiments. When applying forecasting models to imputation tasks, we take the
entire $\mathbf{X}$ as input and force the models to predict the entire
$\tilde{\mathbf{X}}$.

(a) Single Forecasting

(b) Single Imputation

(c) Spatial-Temporal Forecasting

(d) Spatial-Temporal Imputation

Figure 6. Ablation study on the traffic dataset.

###  5.2. Main Results

We compare ReTime with various baselines for single/spatial-temporal time
series forecasting/imputation tasks for different missing rates.

Single Time Series. The results for single time series are presented in Figure
4, where the first and second rows show the Root Mean Squared Error (RMSE) for
forecasting and imputation respectively. For both forecasting and imputation,
compared with RNN based methods, i.e., BRITS and E2GAN, recent block-style
methods, i.e., N-BEATS and Informer, have lower RMSE scores. ReTime has much
lower RMSE scores than N-BEATS and Informer, demonstrating the effectiveness
of the proposed strategy.

An interesting observation for the temperature dataset is that the simple
retrieval baseline REF performs significantly better than other state-of-the-
art baselines, corroborating the power of the reference time series. To be
more specific, firstly, as indicated by the performance of the state-of-the-
art methods, there might exist some complex temporal patterns for the
temperature data, which could not be easily captured by models merely based on
the observed content of the targets. Secondly, the superior performance of REF
over the state-of-the-art methods demonstrates that the retrieved time series
can indeed significantly reduce uncertainty.

(a) Single Forecasting

(b) Spatial-Temporal Forecasting

Figure 7. Effect of using prior snippets for forecasting.

Spatial-Temporal Time Series. The results for spatial-temporal time series are
presented in Figure 5. The general observation is similar to the single time
series that ReTime consistently performs better than the state-of-the-art
methods.

###  5.3. Ablation Study

We study the impact of each component of ReTime on the traffic dataset, which
is the largest dataset.

Impact of the Attentions. As shown in Figure 6, compared with the full model
ReTime , if we only use either the content or temporal attention, the
performance will drop. Besides, spatial attention alone performs worse than
temporal attention alone, showing the importance of modeling temporal
dependencies.

Beyond the First Order Neighbors. In Figure 6, the “1st-order” is ReTime using
the 1st-order neighbors of the targets as references. ReTime is better than
“1st-order”, indicating that it is important to choose appropriate neighbors
as references. RWR used in the relational retrieval stage could capture global
proximity scores.

Effect of the Content Synthesis Model. The “Retrieval Only” in Figure 6 uses
the average of the $K$ retrieved references as the prediction. Compared with
the retrieval-only model, ReTime performs better, demonstrating the
effectiveness of the content synthesis.

Effect of Using Prior Snippets for Forecasting. When the time series has clear
periodic patterns, it is natural to resort to the historical snippets in the
database for help. In Figure 7, “$\Delta T=0$” means that the target and
reference snippets have the same start time, and thus the values of both the
targets and references are zeros after the separation time $\tau$. “$\Delta
T=1$ week” means that the start time of the reference snippets is one week
before the targets. Figure 7 shows that using prior snippets can significantly
improve the model’s performance.

(a) Single Forecasting

(b) Single Imputation

(c) Spatial-Temporal Forecasting

(d) Spatial-Temporal Imputation

Figure 8. The performance of ReTime w.r.t. different $K$.

(a) Time Usage

(b) Memory Usage

Figure 9. Time and Memory Usage.

###  5.4. Impact of $K$

We study the performance of ReTime and its time/memory usage w.r.t. the number
of references $K\in\\{0,1,5,10,20\\}$. The experiments are conducted on the
traffic dataset.

Performance of ReTime. The performance of ReTime w.r.t. the number of
references $K$ is presented in Figure 8, where $r$ denotes the missing rate.
ReTime achieves the best performance when $K\in\\{5,10\\}$ Note that for
spatial-temporal forecasting, the top 1 reference snippet of a target is its
own historical snippet.

Time and Memory Usage. We fix the batch size as 100 and record the average
training time ($\times 10^{-2}$ seconds) of each iteration and the GPU memory
usage ($\times 10^{3}$ MegaBytes) for different $K$. The results in Figure 9
show that the training time and memory usage grow linearly w.r.t. $K$.

##  6\. Related Work

###  6.1. Single Time Series Methods

Many deep learning methods have been proposed to generate time series,
including Recurrent Neural Network (RNN) based methods (Cao et al., 2018a; Che
et al., 2018) and block-style methods (Oreshkin et al., 2020; Zhou et al.,
2021). For example, Cao et al. (Cao et al., 2018a) introduce BRITS which
leverages the recurrent dynamics for both correlated and uncorrelated
multivariate time series. Fortuin et al. (Fortuin et al., 2020) combines the
Gaussian process to capture the temporal dynamics and reconstruct missing
values by VAE (Kingma and Welling, 2013). Luo et al. (Luo et al., 2018)
introduce GRUI and propose a two-stage GAN (Goodfellow et al., 2014) model,
the generator and discriminator of which are based on GRUI. E2GAN (Luo et al.,
2019) further simplifies the generator by combining GRUI and denoising
autoencoder such that the GAN-based imputation can be trained in an end-to-end
manner. Gamboa et al. (Gamboa, 2017) explore different neural networks for
time series analysis. Oreshkin et al. (Oreshkin et al., 2020) introduce
N-BEATS for explainable time series forecasting. Li et al. (Li et al., 2019)
propose an enhanced version of Transformer (Vaswani et al., 2017) for
forecasting. Zhou et al. (Zhou et al., 2021) propose Informer for long time
series forecasting. Wu et al. (Wu et al., 2021) introduce AutoFormer, which
reduces the complexity of Transformer. However, when the missing rate grows to
a high level, the performance of these methods drops significantly. ReTime
address this issue by retrieving relevant reference time series as an
augmentation.

###  6.2. Spatial-Temporal Time Series Methods

Time series often co-evolve with each other. Networks/graphs are commonly used
data structures to model the relations among objects (Tong et al., 2006; Jing
et al., 2021b; Li et al., 2022; Jing et al., 2021c), which have also been used
to model relations in spatial-temporal time series. Traditional methods
leverage probabilistic graphical models to introduce graph regularizations
(Cai et al., 2015a, b). Recently, Li et al. (Li et al., 2018b) introduce DCRNN
combining the diffusion graph kernel with RNN. Yu et al. (Yu et al., 2018) and
Zhao et al. (Zhao et al., 2019) respectively propose STGCN and TGCN for
modeling spatial-temporal time series. Jing et al. (Jing et al., 2021a)
introduce NET3 which captures both explicit and implicit relations among time
series. Cao et al. (Cao et al., 2020) introduce StemGNN which combines graph
and discrete Fourier transform to jointly model spatial and temporal
relations. Wu et al. (Wu et al., 2020) introduce MTGNN which learns the
relation graph of time series. However, these methods do not allow models to
refer to the historical snippets when forecasting future values, which perform
worse than ReTime.

###  6.3. Retrieval Based Generation

The main idea underlying retrieval-based methods is to retrieve references
from databases to guide generation. Cao et al. (Cao et al., 2018b) propose
Re3Sum to generate document summaries based on the retrieved templates. Song
et al. (Song et al., 2018b) propose to generate dialogues based on the
retrieved references. Lewis et al. (Lewis et al., 2020) introduce Retrieval-
Augmented Generation (RAG) for knowledge-intensive natural language processing
tasks e.g., summarization (Jing et al., 2021d). Tseng et al. (Tseng et al.,
2020) propose RetrievalGAN to generate images by retrieving relevant images.
Ordonez et al. (Ordonez et al., 2016) generate image descriptions based on the
retrieved captions. To the best of our knowledge, we present the first
retrieval-based deep generation model for time series data.

##  7\. Conclusion

In this paper, we theoretically quantify the uncertainty of the predicted
values and prove that retrieved references could help to reduce the
uncertainty for prediction results. To empirically demonstrate the
effectiveness of retrieval-based forecasting, we build a simple yet effective
method called ReTime, which is comprised of a relational retrieval stage and a
content synthesis stage. The experimental results on real-world datasets
demonstrate the effectiveness of ReTime.

## References

  * Zhou et al. [2021] Haoyi Zhou, Shanghang Zhang, Jieqi Peng, Shuai Zhang, Jianxin Li, Hui Xiong, and Wancai Zhang.  Informer: Beyond efficient transformer for long sequence time-series forecasting.  In _AAAI_ , 2021.
  * Gamboa [2017] John Cristian Borges Gamboa.  Deep learning for time-series analysis.  _arXiv preprint arXiv:1701.01887_ , 2017.
  * Li et al. [2019] Shiyang Li, Xiaoyong Jin, Yao Xuan, Xiyou Zhou, Wenhu Chen, Yu-Xiang Wang, and Xifeng Yan.  Enhancing the locality and breaking the memory bottleneck of transformer on time series forecasting.  _NeurIPS_ , 2019.
  * Wu et al. [2021] Haixu Wu, Jiehui Xu, Jianmin Wang, and Mingsheng Long.  Autoformer: Decomposition transformers with auto-correlation for long-term series forecasting.  _Advances in Neural Information Processing Systems_ , 34:22419–22430, 2021.
  * Zhou et al. [2020] Dawei Zhou, Lecheng Zheng, Yada Zhu, Jianbo Li, and Jingrui He.  Domain adaptive multi-modality neural attention network for financial forecasting.  In Yennun Huang, Irwin King, Tie-Yan Liu, and Maarten van Steen, editors, _WWW ’20: The Web Conference 2020, Taipei, Taiwan, April 20-24, 2020_ , pages 2230–2240. ACM / IW3C2, 2020.  doi: 10.1145/3366423.3380288.  URL https://doi.org/10.1145/3366423.3380288.
  * Zhu et al. [2022] Yada Zhu, Hanghang Tong, Baoyu Jing, JinJun Xiong, Nitin Gaur, and Anna Wanda Topol.  Network of tensor time series, September 8 2022.  US Patent App. 17/184,880.
  * Yang et al. [2019] Ruihan Yang, Dingsu Wang, Ziyu Wang, Tianyao Chen, Junyan Jiang, and Gus Xia.  Deep music analogy via latent representation disentanglement.  _arXiv preprint arXiv:1906.03626_ , 2019.
  * Li et al. [2018a] Xin Li, Liyu Wu, and Xianfeng Yang.  Exploring the impact of social economic variables on traffic safety performance in hong kong: A time series analysis.  _Safety science_ , 2018a.
  * Sarker et al. [2012] Md Abdur Rashid Sarker, Khorshed Alam, and Jeff Gow.  Exploring the relationship between climate change and rice yield in bangladesh: An analysis of time series data.  _Agricultural Systems_ , 112:11–16, 2012.
  * Wang and Kockelman [2009] Xiaokun Wang and Kara M Kockelman.  Forecasting network data: Spatial interpolation of traffic counts from texas data.  _Transportation research record_ , 2105(1):100–108, 2009.
  * Nalder and Wein [1998] Ian A Nalder and Ross W Wein.  Spatial interpolation of climatic normals: test of a new method in the canadian boreal forest.  _Agricultural and forest meteorology_ , 92(4):211–225, 1998.
  * Bonnet et al. [2001] Philippe Bonnet, Johannes Gehrke, and Praveen Seshadri.  Towards sensor database systems.  In _International Conference on mobile Data management_ , 2001.
  * Pelkonen et al. [2015] Tuomas Pelkonen, Scott Franklin, Justin Teller, Paul Cavallaro, Qi Huang, Justin Meza, and Kaushik Veeraraghavan.  Gorilla: A fast, scalable, in-memory time series database.  _Proceedings of the VLDB Endowment_ , 8(12):1816–1827, 2015.
  * Rhea et al. [2017] Sean Rhea, Eric Wang, Edmund Wong, Ethan Atkins, and Nat Storer.  Littletable: A time-series database and its uses.  In _Proceedings of the 2017 ACM International Conference on Management of Data_ , pages 125–138, 2017.
  * Song et al. [2018a] Dongjin Song, Ning Xia, Wei Cheng, Haifeng Chen, and Dacheng Tao.  Deep r-th root of rank supervised joint binary embedding for multivariate time series retrieval.  In _KDD_ , pages 2229–2238, 2018a.
  * Murphy [2012] Kevin P Murphy.  _Machine learning: a probabilistic perspective_.  MIT press, 2012.
  * Mezard and Montanari [2009] Marc Mezard and Andrea Montanari.  _Information, physics, and computation_.  Oxford University Press, 2009.
  * Tong et al. [2006] Hanghang Tong, Christos Faloutsos, and Jia-Yu Pan.  Fast random walk with restart and its applications.  In _ICDM_ , pages 613–622. IEEE, 2006.
  * Vaswani et al. [2017] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  In _NeurIPS_ , pages 5998–6008, 2017.
  * Ba et al. [2016] Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E Hinton.  Layer normalization.  _arXiv preprint arXiv:1607.06450_ , 2016.
  * Yu et al. [2016] Hsiang-Fu Yu, Nikhil Rao, and Inderjit S Dhillon.  Temporal regularized matrix factorization for high-dimensional time series prediction.  In _NeurIPS_ , 2016.
  * Cao et al. [2020] Defu Cao, Yujing Wang, Juanyong Duan, Ce Zhang, Xia Zhu, Congrui Huang, Yunhai Tong, Bixiong Xu, Jing Bai, Jie Tong, and Qi Zhang.  Spectral temporal graph neural network for multivariate time-series forecasting.  In _NeurIPS_ , 2020.
  * Wu et al. [2020] Zonghan Wu, Shirui Pan, Guodong Long, Jing Jiang, Xiaojun Chang, and Chengqi Zhang.  Connecting the dots: Multivariate time series forecasting with graph neural networks.  In _KDD_ , 2020.
  * Jing et al. [2021a] Baoyu Jing, Hanghang Tong, and Yada Zhu.  Network of tensor time series.  In _The WebConf 2021_ , pages 2425–2437, 2021a.
  * Slivinski et al. [2019] Laura C Slivinski, Gilbert P Compo, Jeffrey S Whitaker, Prashant D Sardeshmukh, Benjamin S Giese, Chesley McColl, Rob Allan, Xungang Yin, Russell Vose, Holly Titchner, et al.  Towards a more reliable historical reanalysis: Improvements for version 3 of the twentieth century reanalysis system.  _Quarterly Journal of the Royal Meteorological Society_ , 145(724):2876–2908, 2019.
  * Oreshkin et al. [2020] Boris N Oreshkin, Dmitri Carpov, Nicolas Chapados, and Yoshua Bengio.  N-beats: Neural basis expansion analysis for interpretable time series forecasting.  _ICLR_ , 2020.
  * Cao et al. [2018a] Wei Cao, Dong Wang, Jian Li, Hao Zhou, Lei Li, and Yitan Li.  Brits: Bidirectional recurrent imputation for time series.  _NeurIPS_ , 2018a.
  * Luo et al. [2019] Yonghong Luo, Ying Zhang, Xiangrui Cai, and Xiaojie Yuan.  E2gan: End-to-end generative adversarial network for multivariate time series imputation.  In _IJCAI_ , pages 3094–3100, 2019.
  * Li et al. [2018b] Yaguang Li, Rose Yu, Cyrus Shahabi, and Yan Liu.  Diffusion convolutional recurrent neural network: Data-driven traffic forecasting.  _ICLR_ , 2018b.
  * Che et al. [2018] Zhengping Che, Sanjay Purushotham, Kyunghyun Cho, David Sontag, and Yan Liu.  Recurrent neural networks for multivariate time series with missing values.  _Scientific reports_ , 8(1):1–12, 2018.
  * Fortuin et al. [2020] Vincent Fortuin, Dmitry Baranchuk, Gunnar Rätsch, and Stephan Mandt.  Gp-vae: Deep probabilistic time series imputation.  In _International Conference on Artificial Intelligence and Statistics_ , pages 1651–1661. PMLR, 2020.
  * Kingma and Welling [2013] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  _arXiv preprint arXiv:1312.6114_ , 2013.
  * Luo et al. [2018] Yonghong Luo, Xiangrui Cai, Ying Zhang, Jun Xu, and Xiaojie Yuan.  Multivariate time series imputation with generative adversarial networks.  In _NeurIPS_ , 2018.
  * Goodfellow et al. [2014] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  In _NeurIPS_ , 2014.
  * Jing et al. [2021b] Baoyu Jing, Chanyoung Park, and Hanghang Tong.  Hdmi: High-order deep multiplex infomax.  In _Proceedings of the Web Conference 2021_ , pages 2414–2424, 2021b.
  * Li et al. [2022] Bolian Li, Baoyu Jing, and Hanghang Tong.  Graph communal contrastive learning.  In _Proceedings of the ACM Web Conference 2022_ , pages 1203–1213, 2022.
  * Jing et al. [2021c] Baoyu Jing, Yuejia Xiang, Xi Chen, Yu Chen, and Hanghang Tong.  Graph-mvp: Multi-view prototypical contrastive learning for multiplex graphs.  _arXiv preprint arXiv:2109.03560_ , 2021c.
  * Cai et al. [2015a] Yongjie Cai, Hanghang Tong, Wei Fan, and Ping Ji.  Fast mining of a network of coevolving time series.  In _SDM_. SIAM, 2015a.
  * Cai et al. [2015b] Yongjie Cai, Hanghang Tong, Wei Fan, Ping Ji, and Qing He.  Facets: Fast comprehensive mining of coevolving high-order time series.  In _KDD_ , 2015b.
  * Yu et al. [2018] Bing Yu, Haoteng Yin, and Zhanxing Zhu.  Spatio-temporal graph convolutional networks: a deep learning framework for traffic forecasting.  In _IJCAI_ , 2018.
  * Zhao et al. [2019] Ling Zhao, Yujiao Song, Chao Zhang, Yu Liu, Pu Wang, Tao Lin, Min Deng, and Haifeng Li.  T-gcn: A temporal graph convolutional network for traffic prediction.  _IEEE Transactions on Intelligent Transportation Systems_ , 21(9), 2019.
  * Cao et al. [2018b] Ziqiang Cao, Wenjie Li, Furu Wei, Sujian Li, et al.  Retrieve, rerank and rewrite: Soft template based neural summarization.  In _ACL_ , 2018b.
  * Song et al. [2018b] Yiping Song, Cheng-Te Li, Jian-Yun Nie, Ming Zhang, Dongyan Zhao, and Rui Yan.  An ensemble of retrieval-based and generation-based human-computer conversation systems.  In _IJCAI_ , 2018b.
  * Lewis et al. [2020] Patrick Lewis, Ethan Perez, Aleksandra Piktus, Fabio Petroni, Vladimir Karpukhin, Naman Goyal, Heinrich Küttler, Mike Lewis, Wen-tau Yih, Tim Rocktäschel, et al.  Retrieval-augmented generation for knowledge-intensive nlp tasks.  In _NeurIPS_ , 2020.
  * Jing et al. [2021d] Baoyu Jing, Zeyu You, Tao Yang, Wei Fan, and Hanghang Tong.  Multiplex graph neural network for extractive text summarization.  In _Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing_ , pages 133–139, 2021d.
  * Tseng et al. [2020] Hung-Yu Tseng, Hsin-Ying Lee, Lu Jiang, Ming-Hsuan Yang, and Weilong Yang.  Retrievegan: Image synthesis via differentiable patch retrieval.  In _ECCV_ , 2020.
  * Ordonez et al. [2016] Vicente Ordonez, Xufeng Han, Polina Kuznetsova, Girish Kulkarni, Margaret Mitchell, Kota Yamaguchi, Karl Stratos, Amit Goyal, Jesse Dodge, Alyssa Mensch, et al.  Large scale retrieval and generation of image descriptions.  In _IJCV_ , 2016.


# Stationary measures on the circle from hyperbolic surfaces with cusps cannot
be straightened by quasi-symmetries

Aitor Azemar  School of Mathematics and Statistics
University of Glasgow
University Place
Glasgow
G12 8QQ United Kingdom<EMAIL_ADDRESS>and  Vaibhav Gadre  School
of Mathematics and Statistics
University of Glasgow
University Place
Glasgow
G12 8QQ United Kingdom<EMAIL_ADDRESS>
###### Abstract.

Stationary measures on the circle that arise from a large class of random
walks on the fundamental group of a finite-area complete hyperbolic surface
with cusps are singular with respect to the Lebesgue measure. In particular,
it is sufficient for singularity that a stationary measure satisfies an
exponential decay for cusp excursions with excursions being measured in the
path metric on horocycles bounding cusps. In this note, we settle a conjecture
of McMullen by proving that the singularity of stationary measures satisfying
such exponential decay is quasi-symmetrically stable, that is under push-
forward by any quasi-symmetry of the circle the measure remains singular.

##  1\. Introduction

Two intervals on the circle are said to be adjacent and equal if they share an
endpoint and have the same length. A circle map is said to be a _quasi-
symmetry_ if there exists $K\geqslant 1$ such that for any pair $I$, $J$ of
adjacent and equal intervals, $\operatorname{Leb}(f(I))\leqslant
K\operatorname{Leb}(f(J))$ and vice versa. Suppose that $P$ is a property that
subsets of the circle might satisfy. We say that a particular set has _quasi-
symmetric stability for $P$_ if its image under any quasi-symmetry satisfies
$P$.

It is of interest to find explicit subsets that display quasi-symmetric
stability. By definition, a quasi-symmetry is Hölder continuous; hence it
preserves zero Hausdorff dimension. However, there are quasi-symmetries that
promote a (Lebesgue) measure zero set to full measure. Specifically, Tukia
defines a quasi-symmetry and a set with Hausdorff dimension less than one such
that its image by the quasi-symmetry has its complement with Hausdorff
dimension less than one. See [21, Main theorem]. By contrast, quasi-symmetries
in higher dimensions are quasi-conformal and hence absolutely continuous.

With the circle being the boundary at infinity for the hyperbolic plane,
singular measures on the circle (and measure zero sets that correspond to
these) arise from various contexts in hyperbolic geometry. A quasi-symmetry of
the circle extends to a quasi-isometry of the hyperbolic plane. Conversely,
quasi-isometries of the hyperbolic plane give rise naturally to quasi-
symmetries of the circle. So it is natural to consider the behaviour of such
measures and sets under quasi-symmetries. In particular, finite area complete
hyperbolic surfaces/orbifolds with cusps have been used to give explicit
measure zero sets that are quasi-symmetrically stable for various properties
(such as measure zero, porosity, etc.). In Section 1.2, we describe some of
these sets and survey known results about them. In this note, we expand the
list by considering singular stationary measures arising from random walks on
_non-uniform lattices_ in $SL(2,\mathbb{R})=\mathrm{Isom}(\mathbb{H}^{2})$
(that is, fundamental groups of finite-area complete hyperbolic
surfaces/orbifolds with cusps).

A random walk on a lattice $\Gamma$ in
$SL(2,\mathbb{R})=\mathrm{Isom}(\mathbb{H}^{2})$ is non-elementary if the
semigroup generated by its support contains a pair of independent hyperbolic
isometries (with distinct stable/unstable fixed points on the circle at
infinity). Projected in to $\mathbb{H}^{2}$ by the group action, a typical
sample path converges to a point on $S^{1}=\partial_{\infty}\mathbb{H}^{2}$.
The convergence defines a stationary measure on
$S^{1}=\partial_{\infty}\mathbb{H}^{2}$ giving the distribution of sample path
limit points.

A lattice $\Gamma$ is non-uniform when the quotient $X=\mathbb{H}^{2}/\Gamma$
is finite area but non-compact and so has cusps. It is well known that
stationary measures from a large class of random walks (for example, with
finite word metric first moment) on a non-uniform lattice are singular. See
[12], [5], [3], [9] for different approaches to understanding the singularity.
A lattice is uniform if the quotient is a closed orbifold. A long standing
conjecture of Guivarc’h–Kaimanovich–Ledrappier states that stationary measures
from finitely supported random walks on a uniform lattice are also singular.
See [5]. Here, we restrict to non-uniform lattices.

For a non-uniform lattice $\Gamma$, we fix cusp neighbourhoods in
$X=\mathbb{H}^{2}/\Gamma$. The lifts of the cusp neighbourhoods give a
horoball packing in $\mathbb{H}^{2}$. We fix a base-point $x$ in
$\mathbb{H}^{2}$. Given a horoball $H$ and a hyperbolic ray $\gamma$ from $x$
that intersects $H$, we define the excursion of $\gamma$ in $H$ to be the
distance (in the path metric on $\partial H$) between the entry and exit
points. Let $r\geqslant 0$ be a constant. We define the $r$-shadow of $H$ to
be the interval at infinity given by hyperbolic rays from $x$ that have an
excursion of size at least $r$ in $H$. We say that a random walk on $\Gamma$
has exponential decay for excursions if for any horoball $H$ and any sample
path converging in to the shadow of $H$, the probability that it converges
into the $r$-shadow decays exponentially in $r$. The decay holds for finitely
supported non-elementary random walks on $\Gamma$. It is expected to be true
for a larger class of random walks (for example with finite word metric first
moment), but this is as yet unproven. We discuss the more general context for
exponential decay in Section 5.4 and sketch a proof for finitely supported
random walks.

Our main result below settles a conjecture of McMullen.

######  Theorem 1.1.

Suppose that $\Gamma$ is a non-uniform lattice in $SL(2,\mathbb{R})$. Suppose
that $\mu$ is a non-elementary probability distribution on $\Gamma$. Let $\nu$
be the stationary measure on $S^{1}$ arising from the $\mu$-random walk.
Suppose that $\nu$ has exponential decay for excursions. Then the pushforward
$f^{*}\nu$ by any quasi-symmetry $f:S^{1}\to S^{1}$ is singular with respect
to the Lebesgue measure on $S^{1}$.

Douady–Earle show that a quasi-symmetry $f:S^{1}\to S^{1}$ always extends to a
bijective quasi-isometry $f$ of $\mathbb{H}^{2}$. See [6, Theorem 2] and also
[13, Theorem 3.1] for a strengthening of it. The pushforward measure
$f^{*}\nu$ is then the stationary measure arising from the same random wlak on
the group of quasi-isometries resulting from conjugating $\Gamma$ by $f$.
Theorem 1.1 is then singularity for stationary measures arising from such
groups of quasi-isometries.

Let $S$ be a surface homeomorphic to $\mathbb{H}^{2}/\Gamma$. We may view the
quotient $\mathbb{H}^{2}/\Gamma$ as defining a marked conformal structure on
$S$ thus giving us a point in the corresponding Teichmüller space
$\mathcal{T}(S)$ of $S$. When $f$ is $\Gamma$-equivariant then the conjugate
$f^{-1}\Gamma f$ is again a non-uniform lattice in $SL(2,\mathbb{R})$
(isomorphic to $\Gamma$). Thus, it can be thought of as a different point of
$\mathcal{T}(S)$. By restricting to $\Gamma$-equivariant quasi-symmetries,
Theorem 1.1 implies the (already known) singularity of stationary measures
given by all of $\mathcal{T}(S)$ from singularity at a single point.

For any finitely supported non-elementary random walk on a uniform lattice
$\Gamma$ in $SL(2,\mathbb{R})$, random walk drift in the hyperbolic plane is a
proper function on the Teichmüller space of the surface that is homeomorphic
to $\mathbb{H}^{2}/\Gamma$. See [2, Theorem A]. It follows that for every such
random walk on the abstract group there is a uniform lattice in
$SL(2,\mathbb{R})$ such that the stationary measure that arises is singular.
If Theorem 1.1 was true for uniform lattices then we would deduce the
singularity conjectured by Guivarc’h–Kaimanovich–Ledrappier. For reasons
outlined in Section 1.5, we find it hard to speculate whether Theorem 1.1
holds for uniform lattices and if it does hold whether proving it is as hard
as proving singularity itself.

###  1.2. Quasi-symmetric stability:

We list some known examples of quasi-symmetric stability.

  * •

Because quasi-symmetries are Hölder continuous, any set with zero Hausdorff
dimension remains so under any quasi-symmetry.

  * •

Any set that is a countable union of uniformly porous sets remains so under
any quasi-symmetry. See [22, Theorem 4.2].

  * •

Any set that is strongly winning or absolutely winning remains so under any
quasi-symmetry. See [17, Theorem 2.2].

###  1.3. Diophantine sets

Some sets arising from geometric considerations have been analysed extensively
for quasi-symmetric stability for properties such as above. The most known
example is the set $D$ of Diophantine directions on a finite area complete
hyperbolic surface $X$ with cusps. $D$ is the set of directions of hyperbolic
rays from a base-point that have bounded excursions in the cusp neighbourhoods
of $X$. By ergodicity of the hyperbolic geodesic flow, it follows that $D$ has
measure zero. Using the concepts of games, McMullen shows that $D$ is an
absolute winning set. McMullen also showed that any absolute winning set
remains so under any quasi-symmetry. See [17, Theorem 1.3 and Theorem 2.2]. It
is known that winning sets have maximal Hausdorff dimension; so $D$ has
Hausdorff dimension one. McMullen also proves that $D$ is countable union of
uniformly porous sets. See [17, Proposition 1.5].

Porosity will be the key concept in our proof of Theorem 1.1. So we give a
proof sketch of the porosity result for $D$. For a positive integer
$m\geqslant 2$, we define $D_{m}\subset D$ consisting of hyperbolic rays whose
excursions are all bounded above by $m$. Given a horoball $H$, the $m$-shadow
of $H$ is a subinterval of the shadow of $H$ with relative size (in the
Lebesgue measure) $\approx 1/m$. The set $D_{m}$ avoids all these sub-
intervals. Using an estimate by Sullivan [20] for the number of horoballs of a
certain size in an interval, one consructs a horoball heirarchy where at each
level we have that

  * •

the horoball shadows have disjoint interiors and cover a set of full measure;
and

  * •

for each horoball of the previous level, the shadows of the current level are
contained in the $m$-shadow or in one of the components (in the horoball
shadow) of the complement of the $m$-shadow.

This shows that the set $D_{m}$ is $1/m$-porous. It follows that $D$ is a
countable union of porous sets, each uniformly porous.

###  1.4. Stationary measures from non-uniform lattices

For the stationary measures, a typical ray has unbounded excursions and a
countable exhaustion by uniformly porous sets is not possible. Using
exponential decay of horoball excursions, we build a horoball heirarchy
$\mathcal{H}^{(n)}$ for which there is a constant $c>0$ (that depends only on
the random walk) such for a typical ray (for the stationary measure) its
excursion in a horoball in $\mathcal{H}^{(n)}$ is eventually (that is, for $n$
large enough depending on the ray) bounded above by $c\log n$. See Proposition
5.7. This implies that eventually a typical ray misses sub-intervals of
horoball shadows at all levels of the heirarchy. The relative sizes of the
sub-intervals in the corresponding shadows approach zero as the level goes to
infinity but quantitatively this happens slowly. This brings us to the notion
of (non-uniform) porosity and we show that the support of the stationary
measure is contained in a countable union of such porous sets.

By a classical theorem of Ahlfors–Beurling (see [1, Theorem 1]), any quasi-
symmetry of the circle is given by integration with respect to a doubling
measure. By a theorem of Wu (see [23, Theorem 1]), non-uniform porous sets
with a high enough porosity (slow enough decay of the relative sizes of
complementary subintervals/holes) are measure zero for any doubling measure,
hence quasi-symmetrically stable. Our excursion upper bounds imply the slow
enough decay of the relative sizes of avoided subintervals required in Wu’s
theorem. Our main theorem, namely Theorem 1.1, is then a direct consequence of
Wu’s theorem.

Singularity of stationary measures is also known in the following analogous
contexts.

  * •

Non-uniform lattices in the isometry group of $\mathbb{H}^{n}$ for $n>2$ (see
[18]);

  * •

Mapping class groups of finite type surfaces with negative Euler
characteristic (see [7], [10]).

The random walk boundaries in the above contexts are the spheres at infinity
$S^{n-1}=\partial_{\infty}\mathbb{H}^{n}$ in the first case, and the Thurston
boundaries of the corresponding Teichmüller spaces (see [15]). The dimensions
of these boundaries are at least two. Since quasi-symmetries (that is, quasi-
conformal maps) in higher dimensions are absolutely continuous, quasi-
symmetric stability holds directly. However, the question of porosity is still
of interest. We expect that the stationary measure has support contained in a
countable union of highly porous (but not uniformly porous) subsets of the
boundaries.

###  1.5. Dimensions of stationary measures:

Our countable collection of sets that have high enough porosity is constructed
using excursion upper bounds. In particular, (up to excising a countable
subset of $S^{1}$) their union contains the set $D$ of Diophantine directions.
Thus, the Hausdorff dimension of their union is one.

On the other hand, the dimension of a singular stationary measure is expected
to be strictly less than one thus giving us a stronger form of singularity.
See [5, Conjecture 1.21] and in a more general form in [14]. For uniform
lattices, finite first moment for a word metric is the same as finite first
moment for the hyperbolic metric. So the dimension gap is conjectured for
finitely supported random walks on uniform lattices, the general picture for
diffused random walks being unclear.

For non-uniform lattices, finite first moment in the word metric implies
finite first moment for the hyperbolic metric but the converse is not true
because the parabolic subgroups (corresponding to the cusps) are exponentially
distorted. Hence, the dimension gap may be conjectured for more general random
walks that have finite word metric first moment. Our countable union of porous
sets while serving the purpose for Theorem 1.1 is nevertheless too large to
give information regarding the dimension of the stationary measure.

Despite the dimension conjecture for finitely supported random walks on
uniform lattices, we think it to be likely that the stationary measure
exhibits no porosity. This should rule out our approach in proving Theorem 1.1
for uniform lattices. It is also possible that there exists a non-equivariant
quasi-symmetry that does straighten the stationary measure in this case.

###  1.6. Bounds on maximal excursions of typical geodesics:

Similar upper bounds hold when all cusp excursions of a typical hyperbolic ray
for the stationary measure are considered. Specifically, there is a constant
$c>0$ that depends only on the random walk such that along a typical ray the
$n^{\text{th}}$ excursion $E_{n}$ is eventually (that is, for all $n$ large
enough depending on the ray) at most $c\log n$. This result follows directly
from the exponential decay of excursions and the easy direction of the
Borel–Cantelli lemma. In fact, the upper bounds also hold in analogous
contexts such as for non-uniform lattices in $\text{Isom}(\mathbb{H}^{n})$ for
$n>2$ and mapping class groups of finite type surfaces with negative Euler
characteristic.

Our construction of the horoball hierarchy is more involved because it is
built to track porosity. So it is not $\Gamma$-equivariant. Excursions in
these horoballs are only a subset of all cusp excursions. Similarity of the
bounds suggest that the set of excursions in the horoballs of the heirarchy
has positive density among all excursions.

The derivation of a logarithmic lower bound for the maximal excursion over all
excursions up to time $T$ is harder and we do not attempt it here. It should
require the following two properties.

  * •

A lower bound for the tail distribution (in the stationary measure) of an
excursion.

  * •

A quasi-independence of excursions.

We expect these properties to be true when the random walk generates $\Gamma$
as a semi-group and is finitely supported (more generally, has a finite first
word metric moment). The arguments are likely to be quite technical. By
knowing logarithmic lower bounds we would describe more precisely a full
measure set for the stationary measure (in particular, it would not contain
$D$). It might be sufficient to prove the dimension gap (referred to in
Section 1.5).

For the Lebesgue measure, the excursions are much larger. In particular, there
is a constant $c^{\prime}>0$ such that for any $c>1$ the largest excursion
$E_{\max}(n)$ along any Lebesgue-typical ray is eventually (that is for $n$
large enough depending on the ray) at most $c^{\prime}n(\log n)^{c}$. See [8,
Lemma 5.5]. Lower bounds for $E_{\max}(n)$ are also known; for example
Sullivan’s $\log$-law says that the $\lim(\log E_{\max}(n)/\log n)$ is $1/\log
2$. See [20].

###  1.7. Acknowledgements

We thank C. McMullen for his comments on a preliminary draft of the paper.

##  2\. Horoball heirarchies

Let $\Gamma$ be a non-uniform lattice in $SL(2,\mathbb{R})$. Let
$N_{1},\cdots,N_{k}$ be (pairwise disjoint) cusp neighbourhoods in the finite
area hyperbolic surface $X=\mathbb{H}/\Gamma$. Lifting these neighbourhoods to
$\mathbb{H}$, let $\mathcal{H}$ be the resulting horoball collection. We call
the complement $\mathbb{H}-\cup_{H\in\mathcal{H}}H$ the _thick part_
$\mathcal{T}$. We fix a base-point $x$ in $\mathcal{T}$.

Let $p,p^{\prime}$ be a pair of distinct points in
$S^{1}=\partial_{\infty}\mathbb{H}$. They partition $S^{1}$ into a pair of
intervals with disjoint interiors. Using the orientation on $S^{1}$
(compatible with orientation on $\mathbb{H}$), we define
$[p,p^{\prime}]\subset S^{1}$ to be the interval that is positively oriented
as we traverse from $p$ to $p^{\prime}$.

Let $H$ be any horoball in $\mathcal{H}$ with $p_{\infty}(H)$ in $S^{1}$ its
point at infinity. Let $\gamma^{-}(H)$ and $\gamma^{+}(H)$ be hyperbolic rays
from $x$ with endpoints $p^{-}(H)$ and $p^{+}(H)$ in $S^{1}$ such that

  * •

$\gamma^{-}(H)$ and $\gamma^{+}(H)$ and tangent to $H$; and

  * •

$p_{\infty}(H)$ is contained in $[p^{-}(H),p^{+}(H)]$.

We call the interval $\operatorname{Sh}_{0}(H)=[p^{-}(H),p^{+}(H)]$ the
_shadow at infinity_ of $H$.

Let $\gamma$ be a geodesic ray from $x$ that converges to a point in
$\operatorname{Sh}_{0}(H)$. Let $x_{p}$ and $x^{\prime}_{p}$ be the points of
intersection of $\gamma$ with $\partial H$. We define the _excursion_
$E(\gamma,H)$ to be the distance along $\partial H$ between $x_{p}$ and
$x^{\prime}_{p}$.

To define shadows more generally, suppose that $r\geqslant 0$. Let
$\gamma^{-}_{r}$ and $\gamma^{+}_{r}$ be geodesic rays from $x$ such that

  * •

$\gamma^{-}_{r}$ converges to a point $p^{-}_{r}$ in
$[p^{-}(H),p_{\infty}(H)]$;

  * •

$\gamma^{+}_{r}$ converges to a point $p^{+}_{r}$ in
$[p_{\infty}(H),p^{+}(H)]$; and

  * •

$E(\gamma^{-}_{r},H)=r$ and $E(\gamma^{+}_{r},H)=r$.

We define the shadow $\operatorname{Sh}_{r}(H)$ to be the interval
$[p^{-}_{r},p^{+}_{r}]$. By definition, if $r<s$ then
$\operatorname{Sh}_{s}(H)$ is _nested_ in $\operatorname{Sh}_{r}(H)$ that is,
the closure of $\operatorname{Sh}_{s}(H)$ is contained in the interior of
$\operatorname{Sh}_{r}(H)$. Let $\operatorname{Leb}$ denote the Lebesgue
measure on $S^{1}$. It follows from the nesting (and basic hyperbolic
geometry) that there is a constant $c>0$ such that for any horoball $H$ in
$\mathcal{H}$

(2.1)
$\operatorname{Leb}(\operatorname{Sh}_{r}(H))\geqslant\frac{c}{r}\operatorname{Leb}(\operatorname{Sh}_{0}(H))$

See [8, Section 3].

We organise horoballs by the (Euclidean) sizes of their shadows as follows.
Suppose that $0<\rho<1$ is a positive constant. We say that a horoball has
$\rho$ size $n$ if
$\rho^{n+1}<\operatorname{Leb}(\operatorname{Sh}_{0}(H))\leqslant\rho^{n}$.
Suppose that $I$ is an interval in $S^{1}$. We define $\mathcal{H}_{n}(I)$ to
be the collection of those horoballs $H$ in $\mathcal{H}$ such that

  * •

$\operatorname{Sh}_{0}(H)$ is contained in $I$; and

  * •

the $\rho$-size of $H$ is $n$.

We recall [20, Proposition 4]. The proof uses the mixing of the hyperbolic
geodesic flow on $T^{1}\mathbb{H}^{2}/\Gamma$. The proposition is not stated
as precisely in [20] as it is stated here. However, our statement is implicit
in the proof of [20, Proposition 4] requiring no modification other than
tracking explicitly how the counting depends on the compact subset of the
boundary.

######  Proposition 2.2.

Suppose that $\Gamma$ is a non-uniform lattice in $SL(2,\mathbb{R})$. There
are positive constants $C$ and $\rho<1$ such that for any interval $I$ in
$S^{1}$ there is a positive integer $n_{I}$ such that

$|\mathcal{H}_{n}(I)|\geqslant C\frac{\operatorname{Leb}(I)}{\rho^{n}}$

for all $n\geqslant n_{I}$.

Suppose that $I$ is an interval in $S^{1}$. We define the subcollection
$\mathcal{H}^{\prime}_{n}(I)\subset\mathcal{H}_{n}(I)$ by requiring that for
any pair of horoballs $H_{1}$ and $H_{2}^{\prime}$ in
$\mathcal{H}^{\prime}_{n}(I)$ the shadows $\operatorname{Sh}_{0}(H_{1})$ and
$\operatorname{Sh}_{0}(H_{2})$ have disjoint interiors. We use Proposition 2.2
to deduce:

######  Lemma 2.3.

Suppose that $\Gamma$ is a non-uniform lattice in $SL(2,\mathbb{R})$. There is
a positive constant $0<\delta<1$ such that for any interval $I$ in $S^{1}$
there is a positive integer $n$ such that the union
$U=\cup_{H\in\mathcal{H}^{\prime}_{n}(I)}\operatorname{Sh}_{0}(H)$ satisfies

$\operatorname{Leb}(U)>\delta\operatorname{Leb}(I)$

###### Proof.

By Proposition 2.2, we may fix $n_{1}=n_{I}$ and consider horoballs $H$ such
that $\operatorname{Sh}_{0}(H)$ is contained in $I$ and the $\rho$-size of $H$
is $n_{1}$. By definition, each horoball shadow can intersect at most $1/\rho$
shadows of size $n$. Thus, we may retain at least
$C^{\prime}\operatorname{Leb}(I)/\rho^{n}$ horoballs (where
$C^{\prime}=c\rho$) to get the collection $\mathcal{H}^{\prime}_{n}(I)$. It
follows that for
$U=\cup_{H\in\mathcal{H}^{\prime}_{n}(I)}\operatorname{Sh}_{0}(H)$ we have

$\operatorname{Leb}(U)=\sum_{H\in\mathcal{H}^{\prime}_{n}(I)}\operatorname{Sh}_{0}(H)\geqslant(|\mathcal{H}^{\prime}_{n}(I)|)\cdot\rho^{n+1}=C\rho^{2}\operatorname{Leb}(I)$

The lemma follows by choosing $\delta<C\rho^{2}$. ∎

We use Lemma 2.3 to deduce:

######  Lemma 2.4.

Suppose that $I$ is an interval in $S^{1}$. There exists a collection
$\mathcal{H}(I)\subset\mathcal{H}$ of horoballs $H$ such that

  * •

for distinct horoballs $H$ and $H^{\prime}$ in $\mathcal{H}(I)$, the interiors
of $\operatorname{Sh}_{0}(H)$ and $\operatorname{Sh}_{0}(H^{\prime})$ are
disjoint; and

  * •

the union $U=\bigcup_{H\in\mathcal{H}(I)}\operatorname{Sh}_{0}(H)$ has full
$\operatorname{Leb}$-measure in $I$.

###### Proof.

By Lemma 2.3 applied to $I$, we can find a collection
$\mathcal{H}(I,1)=\mathcal{H}^{\prime}_{n}(I)$ of horoballs all of $\rho$-size
$n$ such that their shadows have disjoint interiors and their union $U(I,1)$
satisfies $\operatorname{Leb}(U(I,1))>\delta\operatorname{Leb}(I)$. We now
consider the intervals that are the components of $I-U(I,1)$. These are
finitely many and again by Lemma 2.3 applied to each of these intervals, we
can find a collection of horoballs for each interval such that their shadows
have disjoint interiors and their union has proportion (for the Lebesgue
measure) at least $\delta$ in the interval. We may then define
$\mathcal{H}(I,2)$ to be the collection of these horoballs and define $U(I,2)$
to be the union of their shadows. Note that
$\operatorname{Leb}(U(I,2))>\delta\operatorname{Leb}(I-U(I,1))$.

Iterating infinitely, we get a sequence horoball collections
$\mathcal{H}(I,k)$ and subsets $U(I,k)$ of $I$ such that

  * •

$U(I,0)$ is the empty set,

  * •

$U(I,k)\subset I-\cup_{i\leqslant k-1}U(I,i)$;

  * •

for $k\geqslant 1$ we have that $U(I,k)$ is the finite union of shadows of
horoballs in $\mathcal{H}(I,k)$ with pairwise disjoint interiors, and

  * •

there is a constant $\delta>0$ such that
$\operatorname{Leb}(U(I,k))>\delta\operatorname{Leb}(I-\cup_{i\leqslant
k-1}U(I,i))$ for all $k\geqslant 1$.

It follows that the infinite union $\cup_{k}U(I,k)$ is a countable union of
horoball shadows with pairwise disjoint interiors and it has full measure in
$I$. ∎

Using Lemma 2.4, we define the first sub-collection $\mathcal{H}^{(1)}$ of
horoballs in our hierarchy.

######  Lemma 2.5.

There exists a collection $\mathcal{H}^{(1)}\subset\mathcal{H}$ of horoballs
$H$ such that

  * •

for distinct horoballs $H$ and $H^{\prime}$ in $\mathcal{H}^{(1)}$, the
interiors of $\operatorname{Sh}_{0}(H)$ and
$\operatorname{Sh}_{0}(H^{\prime})$ are disjoint; and

  * •

the union $U^{(1)}=\bigcup_{H\in\mathcal{H}^{(1)}}\operatorname{Sh}_{0}(H)$
has full $\operatorname{Leb}$-measure in $S^{1}$.

###### Proof.

The proof follows directly from Lemma 2.4 applied to $S^{1}$. ∎

Let $\alpha>0$ be a constant. The horoball sub-collections of $\mathcal{H}$
that we define below depend on $\alpha$. We drop the dependence from our
notation to keep it simple.

We define the sequence $r_{n}=\alpha\log n$. We now define the collection
$\mathcal{H}^{(n)}$ for $n\geqslant 2$ inductively. Let $H$ be a horoball in
$\mathcal{H}^{(n-1)}$. Let $\mathcal{K}$ be the collection of horoballs
$H^{\prime}$ such that $\operatorname{Sh}_{0}(H^{\prime})$ is contained in
either $\operatorname{Sh}_{r_{n}}(H)$ or in the closure of one of the
components of $\operatorname{Sh}_{0}(H)-\operatorname{Sh}_{r_{n}}(H)$. By
Lemma 2.4, we can find a collection $\mathcal{H}^{(n)}_{H}$ of horoballs in
$\mathcal{K}$ such that

  * •

for distinct horoballs $K$ and $K^{\prime}$ in $\mathcal{H}^{(n)}_{H}$, the
interiors of $\operatorname{Sh}_{0}(K)$ and
$\operatorname{Sh}_{0}(K^{\prime})$ are disjoint; and

  * •

the union $\bigcup_{K\in\mathcal{H}^{(n)}_{H}}\operatorname{Sh}_{0}(K)$ has
full measure in $\operatorname{Sh}_{0}(H)$.

See Figure 2.6.

$p_{\infty}$$p^{-}$$p^{+}$$p^{-}_{r}$$p^{+}_{r}$$\gamma^{-}$$\gamma^{-}_{r}$$\gamma^{+}_{r}$$\gamma^{+}$

Figure 2.6. Illustration of how the collection $\mathcal{H}^{(n)}_{H}$ is
constructed.

We define $\mathcal{L}^{(n)}_{H}$ to be the collection of those horoballs
$H^{\prime}$ in $\mathcal{H}^{(n)}_{H}$ such that
$\operatorname{Sh}_{0}(H^{\prime})$ is contained in
$\operatorname{Sh}_{r_{n}}(H)$. We define $L^{(n)}_{H}$ to be the union
$\bigcup_{H^{\prime}\in\mathcal{L}^{(n)}_{H}}\operatorname{Sh}_{0}(H^{\prime})$.
By Equation 2.1,

$\operatorname{Leb}\left(L^{(n)}_{H}\right)\geqslant\frac{c}{r_{n}}\operatorname{Leb}(\operatorname{Sh}_{0}(H))$

We define

$\mathcal{H}^{(n)}=\bigcup_{H\in\mathcal{H}^{(n-1)}}\mathcal{H}^{(n)}_{H}$

Note that the union
$U^{(n)}=\bigcup_{H\in\mathcal{H}^{(n)}}\operatorname{Sh}(H)$ is full
$\operatorname{Leb}$-measure in $S^{1}$.

We also define

$\mathcal{L}^{(n)}=\bigcup_{H\in\mathcal{H}^{(n-1)}}\mathcal{L}^{(n)}_{H}$

We define $L^{(n)}$ to be the union
$\bigcup_{H\in\mathcal{L}^{(n)}}\operatorname{Sh}_{0}(H)$. Since $U^{(n)}$ has
full measure, by Equation 2.1 we get

$\operatorname{Leb}\left(L^{(n)}\right)\geqslant\frac{c}{r_{n}}\operatorname{Leb}\left(U^{(n)}\right)=\frac{c}{r_{n}}$

##  3\. Porous Sets

We fix $n\geqslant 1$. We define the set $Y(n,n+1)$ by

(3.1) $Y(n,n+1)=U^{(n)}-L^{(n)}$

Let $m>n$. We define the set $Y_{n,m}$ inductively by
$Y(n,m)=Y(n,m-1)-L^{(m-1)}\cap Y(n,m-1)$. We define $X_{n}$ to be the infinite
intersection $\bigcap_{m>n}Y(n,m)$.

######  Definition 3.2.

Let $a_{n}$ be a sequence of numbers satisfying $0<a_{n}<1$. A set $X$ is said
to be $\\{a_{n}\\}$-porous if there is a sequence of coverings
$\mathcal{E}^{(n)}=\\{E(n,j)\\}$ by intervals with disjoint interiors such
that

  * •

$E(n,j)-X$ contains an interval $E^{\prime}(n,j)$ such that
$\operatorname{Leb}(E^{\prime}(n,j)\geqslant a_{n}\operatorname{Leb}(E(n,j))$,

  * •

$\bigcup_{\mathcal{E}^{(n+1)}}E(n+1,k)$ is contained in
$\bigcup_{\mathcal{E}^{(n)}}E(n,j)-E^{\prime}(n,j)$, and

  * •

$\sup_{\mathcal{E}^{(n)}}\operatorname{Leb}(E(n,j))$ tends to $0$ as
$n\to\infty$.

We call the intervals $E^{\prime}(n,j)$ the _holes at level $j$_.

######  Proposition 3.3.

The set $X_{n}$ is porous for the sequence

$a_{j}=\frac{c}{r_{n+j}}$

###### Proof.

For each $j\geqslant 1$, we set $\mathcal{E}^{(j)}$ to be collection of
shadows $\operatorname{Sh}_{0}(H)$ where $H$ is contained in $Y(n,n+j-1)$. For
each such shadow, we set the hole to be $\operatorname{Sh}_{r_{n+j}}(H)$. By
Equation 2.1, the required estimate for the relative measure follows. ∎

##  4\. Quasi-symmetries and Measure zero sets

We say that two intervals $I$ and $J$ in $S^{1}$ are _adjacent and equal_ if
they share an endpoint and have equal length.

######  Definition 4.1.

A map $f:S^{1}\to S^{1}$ is said to be a _quasi-symmetry_ if there exists a
constant $K\geqslant 1$ such that for any adjacent and equal intervals $I$ and
$J$ in $S^{1}$ we have

$\frac{1}{K}\leqslant\frac{\operatorname{Leb}(f(I))}{\operatorname{Leb}(f(J))}\leqslant
K$

######  Definition 4.2.

A measure $\mu$ on $S^{1}$ is said to be a _doubling_ measure if there exists
a constant $M>0$ such that for any $p\in S^{1}$ and any radius $R>0$, the
interval $I(p,2R)$ centred at $p$ and with radius $2R$ satisfies

$0<\mu(I(p,2R))<M\mu(I(p,R))$

By a classical result of Ahlfors–Beurling, namely [1, Theorem 1], a map
$f:S^{1}\to S^{1}$ is quasi-symmetric if and only if it is given by
integration with respect to a doubling measure, that is, if there exists a
doubling measure $\mu$ and $\theta>0$ such that

$f(p^{\prime})-f(p)=\int_{p}^{p^{\prime}}d\mu$

We begin by recalling [23, Theorem 1]

######  Theorem 4.3.

Suppose that $a_{n}$ is a sequence of numbers satisfying $0<a_{n}<1$. Suppose
that for all real numbers $\kappa\geqslant 1$ the sum

$\sum_{n=1}^{\infty}(a_{n})^{\kappa}$

diverges. Then any $\\{a_{n}\\}$-porous set has measure zero for any doubling
measure.

In other words, if the sequence $\\{a_{n}\\}$ satisfies the hypothesis then
for any quasi-symmetric map $f:S^{1}\to S^{1}$ the image of an
$\\{a_{n}\\}$-porous set has measure zero for the Lebesgue measure.

Recall our definition of the set $X_{n}$. By Theorem 4.3, we deduce:

######  Proposition 4.4.

Suppose that $f:S^{1}\to S^{1}$ is a quasi-symmetry. We have

$\operatorname{Leb}(f(X_{n}))=0$

###### Proof.

Note that for any $\kappa\geqslant 1$, the sequence $a_{j}$ in Proposition 3.3
satisfies

$\sum_{j}a_{j}^{\kappa}=\sum_{j}\left(\frac{c}{\alpha\log(n+j)}\right)^{\kappa}=\infty$

The proposition follows from Theorem 4.3. ∎

##  5\. Stationary measures and horoball excursions

###  5.1. Horoball excursions and nested sets in the thick part

The action of the (virtually free) group $\Gamma$ on the thick part
$\mathcal{T}$ is co-compact. By the S̆varc-Milnor lemma, the thick part
$\mathcal{T}$ with the induced path metric $D$ is quasi-isometric to $\Gamma$
for any word metric on $\Gamma$ defined by a finite symmetric generating set.
Let $\pi\colon\mathbb{H}\to\mathcal{T}$ be the closest point projection. We
call the projection $\pi(\gamma)$ of a hyperbolic geodesic $\gamma$ a
_projected path_.

Let $x$ be a point in $\mathcal{T}$ and $H$ a horoball with the point at
infinity $p_{\infty}(H)\in S^{1}$. With $x$ as a base-point, let
$\gamma_{\infty}(H)$ be the geodesic ray from $x$ to $p_{\infty}(H)$. Recall
that $\gamma^{-}(H)$ and $\gamma^{+}(H)$ are rays from $x$ tangent to $H$.
Also recall that $\gamma^{-}_{r}$ and $\gamma^{+}_{r}$ are rays from $x$ that
have excursions of size $r$ in $H$.

Let $s^{-}_{r}$ denote the arc $\pi(\gamma^{-}_{r})\cap\partial H$. Let
$x^{-}_{1}$ and $x^{-}_{2}$ be the endpoints of $s^{-}_{r}$ with $x^{-}_{r}$
the farthest (among the two) from $x$. For any geodesic ray $\gamma$ from $x$
such that $\gamma^{-}_{r},\gamma,\gamma_{\infty}$ are in the counter-clockwise
order, the projected path $\pi(\gamma)$ contains the arc $s^{-}_{r}$. For such
a ray, let $t^{-}_{r}$ denote the initial segment of $\pi(\gamma)$ from $x$ to
$x^{-}_{2}$. Let $R^{-}(r)$ be the set of such rays and define the _thick
shadow_

$C^{-}(x,H;r)=\bigcup_{\gamma\in R^{-}(r)}\pi(\gamma)-t^{-}_{r}$

We similarly define $s^{+}_{r}$ and use it to define the thick shadow
$C^{+}(x,H;r)$ using the ray $\gamma^{+}_{r}$.

We define $\partial_{\infty}C^{-}(x,H;r)$ to be the set of endpoints at
infinity of the rays in $R^{-}(r)$. It follows that
$\partial_{\infty}C^{-}(x,H;r)=[p^{-}_{r}(H),p_{\infty}(H)]$ and
$\partial_{\infty}=[p_{\infty}(H),p^{+}_{r}(H)]$ and so their union is
$\operatorname{Sh}_{r}(H)$.

######  Lemma 5.2.

There exists constants $K\geqslant 1$ and $r_{0}>0$ that depend only on
$\Gamma$ and the choice of cusp neighbourhoods such that for all $r$ large
enough and $(r^{\prime}-r)/K-r_{0}>0$, we have

(5.3)
$D(\mathcal{T}-C^{-}(x,H;r),C^{-}(x,H;r^{\prime}))=\frac{r^{\prime}-r}{K}-r_{0},$

###### Proof.

By [9, Lemma 2.1], a projected path is a quasi-geodesic in $\mathcal{T}$.
Since $C^{-}(x,H;R)$ is defined with a projected arc as a prefix, the lemma
follows. ∎

###  5.4. Exponential decay for cusp excursions

Let $\mu$ be a non-elementary probability distribution on $\Gamma$. We
consider the random walk on $\Gamma$ given by $\mu$. We call the random
product $\omega_{n}=g_{1}\cdots g_{n}$ (where each $g_{i}$ is sampled
independently by $\mu$) a sample path of length $n$. Sample paths of length
$n$ are distributed by the $n$-fold convolution $\mu^{(n)}$. The path space
for the random walk is the space $\Gamma^{\mathbb{N}}$ of one-sided infinite
sequence in $\Gamma$ equipped with the push-forward of the product measure
under the map $(g_{1},g_{2},\cdots)\to(\omega_{1},\omega_{2},\cdots)$. We
denote this measure as $\nu$ and it has marginals $\mu^{(n)}$. The group
$\Gamma$ acts on the left while the random walk increments $g_{n}$ act by
right multiplication $\omega_{n-1}\to\omega_{n}$.

Let $x$ be any base-point in $\mathbb{H}$. For $\nu$-almost every infinite
sample path $\omega^{+}$, the sequence $\omega^{+}_{n}x$ in $\mathbb{H}$
converges to a point $p=p(\omega)$ in $S^{1}$. We may then pushforward $\nu$
to get a measure on $S^{1}$ that gives the distribution of $p$. In a slight
abuse of notation, we continue to denote the measure on $S^{1}$ by $\nu$.

We now define what it means for a random walk to have exponential decay for
excursions.

######  Definition 5.5.

A random walk has _exponential decay for excursions_ if there exists constants
$c\geqslant 0$ and $\alpha>0$ such that for any horoball $H$ in $\mathcal{H}$
and all constants $r,s\geqslant 0$ such that $s\geqslant r+c$, we have

$\nu(\operatorname{Sh}_{s}(H))\leqslant
e^{-\alpha}\nu(\operatorname{Sh}_{r}(H)).$

Note that the exponential decay implies that for $r\geqslant c$ and
$\beta=\alpha((1/c)-(1/r))$

$\nu(\operatorname{Sh}_{r}(H))\leqslant e^{-\beta
r}\nu(\operatorname{Sh}_{0}(H))$

For completeness, we include a proof sketch that finitely supported random
walks that generate the lattice as a semi-group have the desired exponential
decay for excursions. As a special case of more general exponential decay of
shadows (see [16, Lemma 2.10]), the exponential decay is true for any finitely
supported non-elementary random walk. In the specific contexts mentioned here,
we expect the exponential decay to be true for all non-elementary random walks
with finite word metric first moment, but this is open at present.

###### Proof sketch for finitely supported random walks that generate
$\Gamma$:.

For any pair of points $z\neq z^{\prime}$ in $\mathcal{T}$, let
$H(z,z^{\prime})$ be the half-space (in the metric $D$ on $\mathcal{T}$) of
points closer to $z^{\prime}$ than $z$. Let $\partial_{\infty}H(z,z^{\prime})$
be the limit set at infinity of $H(z,z^{\prime})$.

Suppose that $H$ is a horoball. Suppose that $\xi>0$ is a constant. Suppose
that $r^{\prime}=r+3\xi$. Recall that $x\in\mathcal{T}$ is our base-point. We
consider the thick shadows $C^{-}(x,H;r^{\prime})\subset C^{-}(x,H;r)$ and
similarly $C^{+}(x,H;r^{\prime})\subset C^{+}(x,H;r)$. We define the set
$C^{-}(x,H;r+\xi)-C^{-}(x,H;r+2\xi)$ to be the _midpoint set_ for the pair
$C^{-}(x,H;r^{\prime})\subset C^{-}(x,H;r)$. We denote the set by
$M^{-}(\xi)$. Similarly, we define the set
$C^{+}(x,H;r+\xi)-C^{+}(x,H;r+2\xi)$ to be the midpoint set for the pair
$C^{+}(x,H;r^{\prime})\subset C^{+}(x,H;r)$. We denote it by $M^{+}(\xi)$.

By Lemma 5.2, there is $\xi^{\prime}>0$ such that for any
$\xi\geqslant\xi^{\prime}$ there is $u=u(\xi)>0$ such that such that for any
point $x^{\prime}$ in the midpoint set $M^{-}(\xi)$ (similarly $M^{+}(\xi)$)
there is a point $x^{\prime\prime}$ with $D(x^{\prime},x^{\prime\prime})=u$
such that $\partial_{\infty}H(x^{\prime},x^{\prime\prime})$ is contained in
the interior of $[p^{-}_{r}(H),p^{-}_{r^{\prime}}(H)]$ (similarly
$[p^{+}_{r^{\prime}}(H),p^{+}_{r}(H)]$).

Since the random walk generates $\Gamma$, there is a constant $\kappa>0$ such
that for any $z$ such that $D(x,z)=u$, we have
$\nu(\partial_{\infty}H(x,z))>\kappa$.

Let $\delta$ be maximum of $d_{\mathcal{T}}(x_{0},gx_{0})$ over $g$ in support
of $\mu$. Suppose that $H$ is a horoball. Recall the constant $r_{0}$ from
Equation 5.3. Let $\xi=\max\\{\xi^{\prime},K(\delta+r_{0})\\}$. Suppose that
$r^{\prime}>r$ be such that $r^{\prime}=r+3\xi$. Then any sample path that
converges to $\operatorname{Sh}_{r^{\prime}}(H)$ intersects $M^{-}(\xi)$ or
$M^{+}(\xi)$. Conditioning on the point of intersection, we conclude that the
probability for a sample path starting at this point to converge in to
$\operatorname{Sh}_{r}(H)-\operatorname{Sh}_{r^{\prime}}(H)$ is at least
$\kappa$. The claimed decay for $\nu$ now follows. ∎

###  5.6. Upper bounds for excursions in the horoball heirarchy

We now return to the horoball hierarchies $\mathcal{H}^{(n)}$ constructed in
Section 2. Suppose that $p$ lies in $S^{1}$. Suppose that $\gamma$ is the
hyperbolic ray from the base point $x$ to $p$. If $\gamma$ has an excursion in
a horoball $H$ in $\mathcal{H}^{(n)}$ then we denote the excursion by
$E^{(n)}(\gamma)$. By definition, for $n\geqslant 2$, the set $L^{(n)}$ is the
set of points $p$ in $S^{1}$ such that the ray $\gamma$ from $x$ to $p$
satisfies $E^{(n-1)}(\gamma)\geqslant r_{n}$.

We prove:

######  Proposition 5.7.

There exist a constant $\alpha>0$ that depends only on $\mu$ such that for the
horoball collections $\mathcal{H}^{(n)}$ and for $\nu$-almost every $p\in
S^{1}$ there is a positive integer $N=N(\gamma)$ such that with
$r_{n}=\alpha\log n$, we have

$E^{(n-1)}(\gamma)<r_{n}$

for all $n>N$.

###### Proof.

We fix $\alpha>0$ to satisfy $\alpha>1/\beta$. We then consider horoball
heirarchies $\mathcal{H}^{(n)}$ corresponding to the sequence
$r_{n}=\alpha\log n$. By the exponential decay of excursion from Definition
5.5, it follows that the

$\nu\left(L^{(n)}\right)<e^{-\beta r_{n}}=e^{-\beta\alpha\log
n}=\frac{1}{n^{\beta\alpha}}$

It follows that

$\sum_{n\geqslant 2}\nu(L^{(n)})<\infty$

By the Borel–Cantelli lemma (the easy direction), it follows that $\nu(\limsup
L^{(n)})=0$. The conclusion follows. ∎

We deduce:

######  Corollary 5.8.

We have

$\nu\left(\bigcup_{n\geqslant 1}X_{n}\right)=1$

###### Proof.

Suppose that $p$ is $\nu$-typical and suppose that $\gamma$ is the hyperbolic
ray from $x$ to $p$. By Proposition 5.7, there exists a positive number $N$
such that $E^{(n-1)}(\gamma)<r_{n}$ for all $n>N$. Thus, $p$ lies in $X_{N}$
and we get our conclusion. ∎

###### Proof of Theorem 1.1.

Set $X=\cup_{n\geqslant 1}X_{n}$. Suppose that $f:S^{1}\to S^{1}$ is a quasi-
symmetry. Consider $f(X)$. By Proposition 4.4, we have

$\operatorname{Leb}(f(X))\leqslant\sum_{n\geqslant
1}\operatorname{Leb}(f(X_{n}))=0$

On the other hand, by definition of a pushforward measure and Corollary 5.8

$f^{*}\nu(f(X))=\nu(X)=1$

Thus $f^{\ast}\nu$ is singular, as desired. ∎

## References

  * [1] Ahlfors, L. and Beurling, A.  The boundary correspondence under quasi-conformal mappings, Acta. Math. 96 (1956) 125-142.
  * [2] Azemar, A., Gadre, V., Gouëzel, S., Haettel, T., Lessa, P. and Uyanik, C.  Random walk speed is a proper function on Teichmüller space, to appear in J. Mod. Dyn.
  * [3] Blachère, S., Haïssinsky, P. and Mathieu, P.  Harmonic measures versus quasi-conformal measures for hyperbolic groups, Ann. Sci. Éc. Norm. Supér. (4) 44 (2011), no. 4, 683-721.
  * [4] Cooper, D.  Quasi-isometries of hyperbolic spaces are almost isometries. Proc. Amer. Math. Soc 123 (1995), no. 7, 2221-2227.
  * [5] Deroin, B., Kleptsyn, V. and Navas, A.  On the question of ergodicity of minimal group actions on the circle, Mosc. Math. J. 9 (2009), no. 2, 263-303.
  * [6] Douady, A. and Earle, C-J.  Conformally natural extension of homeomorphisms of the circle, Acta Math. 157 (1986), no. 1-2, 23-48.
  * [7] Gadre, V.  Harmonic measures for distributions with finite support on the mapping class group are singular, Duke Math J. 163 (2014), no. 2, 309-368.
  * [8] Gadre, V.  Partial sums of excursions along random geodesics and volume asymptotics for thin parts of moduli spaces of quadratic differentials, J. Eur. Math. Soc. 19 (2017), no. 10, 3053-3089.
  * [9] Gadre, V., Maher, J. and Tiozzo, G.  Word length statistics and Lyapunov exponents for Fuchsian groups with cusps, New York J. Math 21 (2015), 511-531.
  * [10] Gadre, V., Maher, J. and Tiozzo, G.  Word length statistics and singularity of harmonic measure, Comment. Math. Helv 92 (2017), no. 1, 1-36.
  * [11] Gehring, F. and Väisälä, J.  Hausdorff dimension and quasi-conformal mappings, J. London Math. Soc. (2) 6 (1973), 504-512.
  * [12] Guivarc’h, Y. and Le Jan, Y.  Asymptotic winding of geodesic flow on modular surfaces and continued fractions, Ann. Sci. Éc. Norm. Supér. (4) 26 (1993), no. 1, 23-50.
  * [13] Ibragimov, Z.  Quasi-isometric extensions of quasi-symmetric mappings of the real line compatible with composition, Ann. Acad. Sci. Fenn. Math. 35 (2010), no. 1, 221-233.
  * [14] Kaimanovich, V. and LePrince, V.  Matrix random products with singular harmonic measure, Geom. Dedicata 150 (2011), 257-279.
  * [15] Kaimanovich, V and Masur, H.  The Poisson boundary of the mapping class group, Invent. Math. 125 (1996), 221-264.
  * [16] Maher, J.  Exponential decay in the mapping class group, J. Lond. Math. Soc. (2) 86 (2012), no. 2, 366-386.
  * [17] McMullen, C.  Winning sets, quasi-conformal maps and Diophantine approximation, Geom. Funct. Anal. 20 (2010), 726-740.
  * [18] Randecker, A. and Tiozzo, G.  Cusp excursion in hyperbolic manifolds and singularity of harmonic measure, J. Mod. Dyn. 17 (2021), 183-211.
  * [19] Series, C.  The modular surface and continued fractions. J. London Math. Soc. (2) 31 (1985), no. 1, 69-80.
  * [20] Sullivan, D.  Disjoint spheres, approximation by imaginary quadratic numbers and the logarithmic law for geodesics, 149 (1982), no. 3-4, 215-237.
  * [21] Tukia, P.  Hausdorff dimension and quasi-symmetric mappings, Math. Scand. 65 (1989), no. 1, 152-160.
  * [22] Väisälä, J.  Porous sets and quasi-symmetric maps, Trans. Amer. Math. Soc. 299 (1987), no. 2, 525-533.
  * [23] Wu, J-M.  Null sets for doubling and dyadic doubling measures, Ann. Acad. Sci. Fenn. Ser. A I Math. 18 (1993), no. 1, 77-91.