# Proximity Graph Maintenance for Fast Online Nearest Neighbor Search

Zhaozhuo Xu, Weijie Zhao, Shulong Tan, Zhixin Zhou, Ping Li

Cognitive Computing Lab
Baidu Research
10900 NE 8th St. Bellevue, WA 98004, USA

{zhaozhuoxu, zhaoweijie12, laos1984, zhixin0825<EMAIL_ADDRESS>
###### Abstract

111Currently, Zhaozhuo Xu is a PhD student at Rice University, Weijie Zhao is
an Assistant Professor at Rochester Institute of Technology, and Zhixin Zhou
is an Assistant Professor at the City University of Hong Kong. The work was
conducted in 2018 and 2019 while Zhaozhuo Xu was a research intern (and a
graduate student at Stanford) and both Zhixin Zhou and Weijie Zhao were
Postdoctoral Researchers, in Cognitive Computing Lab, Baidu Research.

Approximate Nearest Neighbor (ANN) search is a fundamental technique for
(e.g.,) the deployment of recommender systems. Recent studies bring proximity
graph-based methods into practitioners’ attention—proximity graph-based
methods outperform other solutions such as quantization, hashing, and tree-
based ANN algorithm families. In current recommendation systems, data point
insertions, deletions, and queries are streamed into the system in an online
fashion as users and items change dynamically. As proximity graphs are
constructed incrementally by inserting data points as new vertices into the
graph, online insertions and queries are well-supported in proximity graph.
However, a data point deletion incurs removing a vertex from the proximity
graph index, while no proper graph index updating mechanisms are discussed in
previous studies. To tackle the challenge, we propose an incremental proximity
graph maintenance (IPGM) algorithm for online ANN. IPGM supports both vertex
deletion and insertion on proximity graphs. Given a vertex deletion request,
we thoroughly investigate solutions to update the connections of the vertex.
The proposed updating scheme eliminates the performance drop in online ANN
methods on proximity graphs, making the algorithm suitable for practical
systems.

##  1 Introduction

Approximate Nearest Neighbor (ANN) search is a fundamental problem with wide
applications in recommendation systems, search engine, and advertising
(Pazzani and Billsus, 2007; Vanderkam et al., 2013; Xu et al., 2018b; Fan et
al., 2019; Tan et al., 2020b, 2021b; Yu et al., 2022). Given a query vector
$q$ and a dataset $\mathcal{D}$, the goal of ANN search is to find a vector
$p\in\mathcal{D}$ that minimizes metric distance $f(p,q)$ with sub-linear time
complexity. Here, $f$ represents a distance measure, e.g., Euclidean distance,
cosine, Jaccard, GMM (generalized min-max) (Li, 2018) etc. The research on ANN
dated back to the early days of computer science, i.e., KD-tree algorithms
(Friedman et al., 1975, 1977). Even recently, tree-based methods are still
active research topics (Lin et al., 1994; Berchtold et al., 1996; Zezula et
al., 1998; Sakurai et al., 2000; Ram and Sinha, 2019). Another popular line of
research on ANN is based on hashing, including random sampling, (Gaussian and
non-Gaussian) random projections, minwise hashing, consistent weighted
sampling, extremal process, etc. (Broder, 1997; Broder et al., 1998; Indyk and
Motwani, 1998; Gionis et al., 1999; Charikar, 2002; Datar et al., 2004; Li and
Church, 2005; Weiss et al., 2008; Shrivastava and Li, 2012; Wei et al., 2014;
Liong et al., 2015; Gao et al., 2015; Shen et al., 2015; Zhu et al., 2016; Li
and Zhang, 2017; Li, 2019; Li and Li, 2019; Li et al., 2021). Quantization-
based ANN methods are also commonly used (Jégou et al., 2011a; Ge et al.,
2013; Zhang et al., 2014; André et al., 2015; Cao et al., 2016).

In this paper, we focus on graph-based ANN methods (Paredes and Chávez, 2005;
Hajebi et al., 2011; Malkov et al., 2014; Morozov and Babenko, 2018; Fu et
al., 2019; Tan et al., 2019; Zhou et al., 2019; Malkov and Yashunin, 2020;
Zhao et al., 2020; Tan et al., 2021a), which have successfully demonstrates
their impressive empirical performance. Graph-based ANN search methods share
the spirits with idea of six degrees of separation (Karinthy, 1929; Guare,
1990; Kleinberg, 2000). The general paradigm of these methods can be
summarized as two steps: First, we index the dataset as a proximity graph $G$,
where every vertex represents a data vector and every edge connects two
“close” vertices. Then, we walk on $G$ towards the query $q$ to obtain its
nearest neighbors: we maintain a priority queue that stores the searching
vertex candidates. In each search step, we extract the vertex that is closest
to $q$ from the priority queue and push all unvisited neighbors of the vertex
into the priority queue for the future search steps. The proximity graph
indexing is constructed in a bootstrap manner. We insert data vectors as new
vertices incrementally. Each time a new vertex is inserted, we query on the
existing graph index to obtain its nearest neighbors from the priority queue
and create edges between them.

However, graph-based methods do not support ANN search in online settings due
to the lack of vertices update mechanism. In real-world applications, the
online ANN search means that the dataset $D$ is updated by deletion and
insertion for a period of time. For instance, in sponsored search systems for
online advertising, the deep Click-Through Rate (CTR) prediction models (Dave
and Varma, 2010; Fan et al., 2019; Yu et al., 2022) learn user vectors and Ad
vectors and measure them by cosine distance to predict the final CTR. Given a
user vector, ANN search is applied to retrieve ad vectors that minimize their
cosine similarity. Due to the dynamic changes of sponsored items over time,
the online Ads systems will update the ad vectors by deleting the expired Ads
and inserting the recently arrived Ads. In these settings, although the
incremental vertex insertion construction algorithm of graph-based algorithm
can handle the upcoming Ads vectors, there is no specialized design for
removing expired vertices on graph. If we preserve the vertices of expired Ads
on the graph, the search precision and efficiency will drop because a branch
of time is wasted in visiting and ranking useless vertices. If we remove these
vertices as well as their edges, the connectivity of the graph will be broken,
which will also cause inefficiency in the search phase. There are few
literatures on the graph-based online ANN search approaches. Therefore, when
graph-based methods are deployed in online systems, a periodical
reconstruction of indexing graph is inevitable. As re-indexing the large scale
datasets scales as $\mathcal{O}(n\log(n))$ (Malkov et al., 2014), which
increases the computation and memory cost for online services.

To address those challenges, this paper studies online ANN search with graph-
based indexing methods. We imitate the real-world settings by incrementally
performing insertion and deletion in given ANN search vector sets. To
dynamically handle the changes in data distribution, we propose an incremental
proximity graph maintenance (IPGM) algorithm that supports online insertion
and deletion of vertices on graph, which enables the indices update when the
distribution of data is evolving periodically.

Summary of main contributions: Firstly, we formally define the proximity graph
maintenance problem and theoretically prove the feasibility of incremental
updates on proximity graphs. Then, we propose an online graph maintenance
framework and 4 novel proximity graph update strategies that maintain
proximity graph properties when deleting vertices. Massive experimental
evaluations on public datasets confirm the effectiveness—our recommend global
reconnection algorithm outperform the reconstruction baseline in query
processing speed by a maximum factor of 18.8$\%$.

##  2 Related Work

###  2.1 Graph-based ANN Search

Graph-based ANN search methods (Paredes and Chávez, 2005; Hajebi et al., 2011;
Malkov et al., 2014; Morozov and Babenko, 2018; Fu et al., 2019; Tan et al.,
2019; Zhou et al., 2019; Malkov and Yashunin, 2020; Zhao et al., 2020; Tan et
al., 2021a) have received considerable research attention. Most ANN search
scenarios studied by these papers are static. In the experiments section, the
graph-based methods are compared with other ANN search methods on a fixed
dataset such as Sift (Jégou et al., 2011a) or Glove (Pennington et al., 2014).
In recommendations systems, one of the core tasks is: given a user embedding,
to find relevant item embeddings. This is usually realized by learning an
auxiliary data structure concurrently with learning item embeddings (Zhu et
al., 2018, 2019; Zhuo et al., 2020; Gao et al., 2020; Yu et al., 2022). (Tan
et al., 2020a, 2021b) discover that graph-based ANN techniques can be
leveraged on complex neural network similarity measure, i.e., user-item
matching score, to efficiently retrieve recommendation candidates. However, in
current ANN applications like feed-in Ads recommendation (Fan et al., 2019),
the distribution of dataset (e.g., product vectors) changes over time due to
the activation and expiration of products. Few works addressed this problem
except Singh et al. (2021). Therefore, modification to graph-based methods for
online manner is essential for filling their gap with practical usage.
However, as the insertion or deletion of vertex on graph breaks its topology,
direct modification of graph-based methods for online manner is challenging
due to the unaffordable computational cost.

###  2.2 Online ANN Search

Besides graph-based ANN methods, some ANN methodologies can be modified for
online settings, such as hashing methods and quantization methods. However,
their shortcomings lie in the lower speed during the search phase.

Ghashami and Abdullah (2015); Leng et al. (2015); Cakir et al. (2017) modify
the supervised hashing algorithms to accommodate the sequential inputs. These
online hashing models achieve promising results in handling sequential input
vectors. However, these online hashing requires expensive and infeasible
labeling for each data vectors in real-world systems. Meanwhile, these online
learning paradigms do not support the expired data vectors removal.

For quantization methods, the main challenge is to update the codebook
according to the change of data distribution in database. In online settings,
the changes in the distribution in database require expensive re-computation
of codebook and codewords. To tackle this issue, Xu et al. (2018a) propose a
PQ approach that supports data insertion and deletion by updating the codebook
without changing the codewords of existing vectors in the database. This
online PQ method achieves promising results in sliding window settings.

##  3 Problem Statement

Formally, the ANN search problem is defined as follows:

######  Problem 1 (ANN Search).

given a dataset $D$ and a query vector $q$ in Euclidean space $\mathbb{R}^{d}$
we aim at efficiently computing

$\displaystyle p=\arg\max_{x\in D}f(x,q).$ (1)

Here, we discuss metric similarity such as $f(x,q)=-\|x-q\|$, which minimizes
the Euclidean distance, or $f(x,q)=\frac{x^{T}q}{\|x\|\|q\|}$, which maximizes
the cosine similarity. In most real-world datasets, linear scan is infeasible
to solve this problem because $|D|\gg d$. To remedy this issue, we index $D$
as a proximity graph $G=(D,E)$, where each vertex represents a $x_{i}\in D$
and each pair $(x_{i},x_{j})$ indicates the connection between vertices. Given
each $q$, we find the solution of Eq. (1) by walking on $G$. As a walking path
only covers a fraction of vertices on $G$, the speed for solving Eq. (1) is
accelerated at the expense of precision.

For metric similarities that satisfy triangle inequality, an efficient
proximity graph $G$ usually approximates _Delaunay graph_ , which is defined
by _Voronoi cells_.

######  Definition 1.

For fixed $x_{i}\in D\subset\mathbb{R}^{d}$ and a given function $f$, the
Voronoi cells $R_{i}$ is defined as

$\displaystyle R_{i}:=R_{i}(f,D):=\\{q\in\mathbb{R}^{d}\ |\ \forall x\in
D,f(x_{i},q)\geq f(x,q)\\}.$

_Voronoi cells_ determine the solution of Eq. (1). We observe that,
$x_{j}\in\arg\max_{x_{i}\in D}f(x_{i},q)$ if and only if $q\in R_{j}$. The
_Delaunay graph_ is the dual diagram of _Voronoi cells_ , which is given by
the following definition.

######  Definition 2.

For fixed function $f$ and dataset $D\subset\mathbb{R}^{d}$, and given
_Voronoi cells_ $R_{i}$, $i=1,2,\dots,n$ w.r.t. $f$ and $D$, the Delaunay
graph of $D$ is an undirected graph $G$ with vertices $\\{x_{i}\in D\\}$, and
the edge $\\{x_{i},x_{j}\\}$ exists if $R_{i}\cap R_{j}\neq\emptyset$.

For graph-based ANN search problem, the guarantee of finding the optimal
solution for Eq. (1) on _Delaunay graph_ is given in Theorem 1, with proof
provided in Morozov and Babenko (2018).

######  Theorem 1.

For given metric similarity $f$, we assume for any dataset $D$, each the
Voronoi cell $R_{i}$ is a connected. Let $G=(D,E)$ be the _Delaunay graph_
w.r.t. the Voronoi cells. Then for any $q\in\mathbb{R}^{d}$, greedy search on
_Delaunay graph_ returns the solution of (1). In other words, let
$N(x_{i})=\\{x_{j}\in D:\\{x_{i},x_{j}\\}\in E\\}$ be the neighbors of $x_{i}$
on $G$. If $x_{i}$ satisfies

$\displaystyle f(x_{i},q)\geq\max_{x_{j}\in N(x_{i})}f(x_{j},q),$ (2)

then $x_{i}$ is a solution of Eq. (1).

For large scale dataset, the computation cost of constructing and storing both
_Voronoi cells_ and _Delaunay graph_ are unaffordable. The performance of
proximity graphs in ANN search depends on their approximation to _Delaunay
graph_.

In this work, we discuss the maintenance of $G$ in online settings, where the
ANN search problem becomes:

######  Problem 2 (Online ANN Search).

Given a dataset sequence $\\{D_{1},D_{2},\dots D_{T}|D_{i}\cap
D_{i+1}\neq\emptyset,i=1,2,\dots,T\\}$, and a query vector $q$ in Euclidean
space $\mathbb{R}^{d}$ we aim at efficiently computing

$\displaystyle p=\arg\max_{x\in D_{i}}f(x,q)\quad i\in\\{1,2,\dots,T\\}$ (3)

Here, $D_{i}\cap D_{i+1}\neq\emptyset$ indicates the possibility of data
insertion and deletion in each step. Start from proximity graph $G_{1}$ that
approximates _Delaunay graph_ of $D_{1}$, our goal is to develop a proximity
graph maintenance algorithm $\emph{update}(G,D)$ so that
$G_{i+1}=\emph{update}(G_{i},D_{i+1})\quad i\in\\{1,2,\dots,T\\}$ is an
approximation of _Delaunay graph_ with respect to $D_{i+1}$.

##  4 Incremental Proximity Graph

We first theoretically prove that deletion a vertex on Delaunay Graph only
requires update of the removed vertex’s neighbors’ connections. Based on this
insight, we propose a high-level incremental proximity graph maintenance
framework for online update settings.

###  4.1 Theoretical Analysis

In this section, we prove that we only require to update the edges of the
deleted vertex’s neighbors to maintain the properties of Delaunay graph—the
incremental graph maintenance is feasible in online settings. Given the
Delaunay graph $G$ of vertices $S$, suppose we remove an vertex $x_{0}$ from
$S$, and

  1. 1.

suppose $x\in D$ is not a neighbor of $x_{0}$ in $G$, then the neighbors of
$x$ in the new Delaunay graph (corresponding to $S\backslash\\{x_{0}\\}$)
remains the same as the ones in the old $G$.

  2. 2.

suppose $x\in D$ is a neighbor of $x_{0}$ in $G$, then it’s incident edges
except $\\{x,x_{0}\\}$ remains existing in the Delaunay graph.

These properties are rigorously stated in the following theorem. Given a graph
$G$, $G$ is a _Delaunay graph_ or contains _Delaunay graph_ as subgraph, we
demonstrate that

  1. 1.

If a vertex $x_{i}$ is removed from $G$, $G$ still contains a _Delaunay graph_
if we update the connections of $N(x_{i})$.

  2. 2.

If a vertex $x_{i}$ is inserted to $G$, $G$ still contains a _Delaunay graph_
if we set up $N(x_{i})$.

######  Theorem 2.

Let $f$ be continuous function such that for any datasets, every Voronoi cell
is nonempty and connected. Let $G$ and $G^{\prime}$ be the Delaunay graphs
corresponding to data points $D=\\{x_{0},\dots,x_{n}\\}$ and
$D^{\prime}=\\{x_{1},\dots,x_{n}\\}=D\backslash\\{x_{0}\\}$. Let $N(x_{i},G)$
and $N(x_{i},G^{\prime})$ be the set of neighbors of $x_{i}$ on graph $G$ and
$G^{\prime}$ respectively, then

  * (a)

for all $x_{i}\notin N(x_{0},G)$, $N(x_{i},G)=N(x_{i},G^{\prime})$;

  * (b)

removing $x_{0}$ and its the incident edges from $G$, the resulting graph is a
subgraph of $G^{\prime}$.

###### Proof.

(a) We consider two cases.

Firstly, if $N(x_{i},G)\subset N(x_{i},G^{\prime})$, we show that, by the
definition of Voronoi cells, for all $x\in D^{\prime}$ where
$D^{\prime}\subset D$, we have

$\displaystyle R_{x}(D)=$ $\displaystyle\bigcap_{y\in D}\\{q\in Y:f(x,q)\geq
f(y,q)\\}=R_{x}(D^{\prime})$ $\displaystyle\cap\\{q\in Y:f(x,q)\geq
f(x_{0},q)\\}.$

Hence $R_{x}(D)\subset R_{x}(D^{\prime})$. As a result, $R_{x}(D)\cap
R_{y}(D)\neq\emptyset$ implies $R_{x}(D^{\prime})\cap
R_{y}(D^{\prime})\neq\emptyset$. Thus $N(x_{i},G)\subset N(x_{i},G^{\prime})$.

Secondly, if $N(x_{i},G^{\prime})\subset N(x_{i},G)$, suppose $y\in
N(x_{i},G^{\prime})$ but $y\notin N(x_{i},G)$, then by definition of Delaunay
graph, we have

$\displaystyle R_{y}(D^{\prime})\cap
R_{x_{i}}(D^{\prime})\neq\emptyset\quad\text{ and }\quad R_{y}(D)\cap
R_{x_{i}}(D)=\emptyset.$

Since

$\displaystyle R_{y}(D)\cap R_{x_{i}}(D)=$ $\displaystyle
R_{y}(D^{\prime})\cap R_{x_{i}}(D^{\prime})\cap\\{q\in Y:f(y,q)\geq
f(x_{0},q)\\}$ $\displaystyle\cap\\{q\in Y:f(x_{i},q)\geq
f(x_{0},q)\\}=\emptyset$

we have

$\displaystyle R_{y}(D^{\prime})\cap R_{x_{i}}(D^{\prime})\subset$
$\displaystyle\\{\\{q\in Y:f(y,q)<f(x_{0},q)\\}$ $\displaystyle\cup\\{q\in
Y:f(x_{i},q)<f(x_{0},q)\\}\\}$

Hence $R_{y}(D^{\prime})\cap R_{x_{i}}(D^{\prime})\subset R_{x_{0}}(D)$. Pick
$q\in R_{x_{i}}(D)$ and $q^{\prime}\in R_{y}(D^{\prime})\cap
R_{x_{i}}(D^{\prime})$. Since $R_{x_{i}}(D^{\prime})$ is connected, there
exists a path $c:[0,1]\to R_{x_{i}}(D^{\prime})$ connecting $q$ and
$q^{\prime}$. Considering the function

$\displaystyle g(t)=f(x_{0},c(t))-f(x_{i},c(t))$

we have

$\displaystyle g(0)=f(x_{0},q)-f(x_{i},q)\leq 0\qquad
g(1)=f(x_{0},q^{\prime})-f(x_{i},q^{\prime})\geq 0$

By continuity of $f$ and intermediate value theorem, there exists
$t^{*}\in[0,1]$ such that $g(t^{*})=0$. Hence $R_{x_{i}}(D)$ and
$R_{x_{0}}(D)$ intersect at $c(t^{*})$, which contradicts with $x_{i}\notin
N(x_{0},G)$.

(b) Suppose $x,y\in D$ and the edge $\\{x,y\\}\in G$, then $R_{x}(D)\cap
R_{y}(D)\neq\emptyset$. Hence, $R_{x}(D)\cap R_{y}(D)\subset
R_{x}(D^{\prime})\cap R_{y}(D^{\prime})\neq\emptyset$. Therefore,
$\\{x,y\\}\in G^{\prime}$. ∎

It was demonstrated in Navarro (2002) that greed search on proximity graph is
guaranteed to find the exact nearest neighbor for every given query. However,
it is computationally infeasible to build exact Delaunay Graph in high-
dimensional data because the number of edges grows exponentially as the scale
of dataset increases (Beaumont et al., 2007). A typical solution towards this
issue is to approximate Delaunay Graph and trade accuracy for significant
efficiency improvements. Practical applications employ the proximity graph as
an approximation of the Delaunay graph. Therefore, we focus on the incremental
maintenance on proximity graphs in this paper.

###  4.2 High-Level Proximity Graph Maintenance

In this section, we introduce our high-level incremental proximity graph
maintenance algorithm (Algorithm 3).

In each state of the workload, we approximate the Delaunay Graph via a
directed proximity graph, where each vertex have edges pointed to its
neighbors, namely out-neighbors. We use
$W=\\{(o_{1},x_{1}),(o_{2},x_{2}),\dots,(o_{n},x_{n})\\}$ to represent the
online ANN search workload. For each $(o_{i},x_{i})\in W$, $o_{i}$ is the type
of the operations–namely query, insert, or delete–and $x_{i}$ is the
corresponding vectors.

Query. Querying operation is performed by executing the greedy search on the
proximity graph $G$. The algorithm for greedy search is present in Algorithm
1. Given a query vector $q$ and a graph $G$, the algorithm first chooses a
vertex $x_{s}$ randomly on $G$. Starting from $x_{s}$, $q$ walks on $G$.
During the walking, we preserve a priority queue where the vertex that
maximizes the measurement function $f$ is on top of the queue. For every
visited vertex $x$, we push its out-neighbors in the queue after calculating
their similarity with $q$.

1: Input: graph $G=(D,E)$, query vector $q$, candidate size $k$, removed
vertices $Y$.

2: Initialize the priority queue by random sampling,
$S\leftarrow\\{x_{i}|x_{i}\in G,x_{i}\notin Y\\}$.

3: Mark elements in $\\{x_{i}|x_{i}\in G,x_{i}\notin Y\\}$ as unvisited.

4: Mark $x_{i}$ as visited.

5: if $|S|>k$ then

6: $S\leftarrow$ top-$k$ vectors in $S$ that has the largest value of
$-\|x-q\|$. We place them in a descending order.

7: end if

8: while $\exists x\in D$ unvisited and $S$ do

9: $S\leftarrow S\cup\\{y\in D:x\in S,y\text{ unvisited},(x,y)\in E,y\notin
Y\\}$

10: Mark $S$ as visited.

11: if $|S|>k$ then

12: $S\leftarrow$ top-$k$ vectors in $S$ that has the largest value of
$-\|x-q\|$. We place them in a descending order.

13: end if

14: end while

15: Output: $C$.

Algorithm 1 GREEDY-SEARCH$(G,q,k,Y)$

1: Input: vector $x$, $k$-neighbor set $S$ of $x$, out-neighbor degree
threshold $d$, invalid set $I$.

2: Set $N_{x}$ as $\emptyset$.

3: Place $y_{i}\in S$ in descending order of $-\|x-y_{i}\|$.

4: $i\leftarrow 1$.

5: while $|N_{x}|\leq d$ and $i\leq|S|$ do

6: if $\|x-y_{i}\|\leq\min_{z\in N_{x}}\|z-y_{i}\|$ and $y_{i}\notin I$ then

7: $N_{x}\leftarrow N_{x}\cup\\{y_{i}\\}$.

8: end if

9: $i\leftarrow i+1$.

10: end while

11: Output: a set of elements $N_{x}$.

Algorithm 2 SELECT-NEIGHBORS$(x,S,d,I)$

1: Input: workload $W=\\{(o_{1},x_{1}),\cdots,(o_{n},x_{n})\\}$, priority
queue length $k$, out-neighbor degree threshold $d$,measurement function $f$

2: $Y\leftarrow\\{\\}$.

3: for $i=1$ to $n$ do

4: if $o_{i}=query$ then

5: $C_{i}\leftarrow\text{\small GREEDY-SEARCH}(x_{i},G,k,f,Y)$.

6: else if $o_{i}=Insert$ then

7: $C\leftarrow\text{\small GREEDY-SEARCH}(x_{i},G,k,f,Y)$.

8: $N\leftarrow\text{\small SELECT-NEIGHBORS}(x_{i},C,d,\emptyset)$.

9: for $z\in N$ do

10: Add edges $(x_{i},z)$ to $G$ and $(z,x_{i})$ to $G^{\prime}$.

11: end for

12: else if $o_{i}=delete$ then

13: if $x_{i}\in G$ then

14: DELETE-UPDATE-EDGES

115: end if

16: end if

17: end for

18: Output: top-$K$ objects $C_{i}\subset D$ in descending order of $f$ with
$x_{i}$ for every $\\{(o_{i},x_{i})\in W|o_{i}=query\\}$.

Algorithm 3 GRAPH-MAINTENANCE$(W,k,d,f)$

Insertion. When $o_{i}$ is inserted, we perform a greedy search on $G$ to
obtain $x_{i}$’s Top-K neighbors (same as in the query phase). We then extract
the vertex from the queue as the start point for next search step. When the
walking stops, we select the Top-K vertices in the queue to return. The next
step is to select the Top-K neighbors of $x_{i}$ found on $G$ and set
connections in $G$ and $G^{\prime}$. Here we adapt the edge selection
algorithm (Malkov et al., 2014) (Algorithm 2). After the selection, $x_{i}$’s
outgoing vertices are its nearest neighbors in diverse directions.

Deletion. If $o_{i}$ is delete and $x_{i}\in G$, we start the vertex deletion
and graph update process. For every expired vertex $x_{i}$, we first obtain
its incoming vertices on $G$, by the prepared reversed graph $G^{\prime}$. Let
$N^{\prime}(x_{i})$ be the in-neighbors of vertex $x_{i}$ (the vertices that
points to $x_{i}$ in $G$). DELETE-UPDATE-EDGES updates the edges of
$N^{\prime}(x_{i})$ to maintain the proximity graph—the nearest neighbors can
be retrieved from the updated graph by the greedy searching algorithm.
Besides, the reverse graph $G^{\prime}$ is updated according to $G$.

A sub-optimal updating algorithm reduces the connectivity of the proximity
graph and hence the greedy searching cannot locate accurate nearest neighbors.
In the online settings, the deletion operation is called regularly in the
workload—the disadvantage of sub-optimal updating algorithms are magnified.
Therefore, we delve into the design and implementation of the DELETE-UPDATE-
EDGES algorithm in the following paper.

##  5 DELETE-UPDATE-EDGES Algorithms

We investigate online proximity graph deletion algorithms in this section. We
begin with straightforward solutions and then keep refining the update
algorithm step by step.

Figure 1: Illustrations of graph updating algorithms.

###  5.1 Pure Delete

The most straightforward delete operations for proximity graph update is to
remove the vertex as well as its connections. Here we demonstrate our pure
deletion algorithm in Algorithm 4 and Figure 1 (b). Given a deleted vertex
$x_{i}$, we first obtain its out-neighbors and in-neighbors. For each out-
neighbor $x_{j}$, we remove edge $(x_{i},x_{j})$ on proximity graph $G$ and
$(x_{j},x_{i})$ on the reverse proximity graph $G^{\prime}$. For each in-
neighbor $x_{k}$, we remove edge $(x_{k},x_{i})$ on proximity graph $G$ and
$(x_{i},x_{k})$ on the reverse proximity graph $G^{\prime}$. Pure deletion
breaks the connectivity of the proximity graph, which may lead to the fail of
greedy search as the path towards the exact nearest neighbors are blocked.

###  5.2 Vertex Masking

To tackle the preserve the connectivity challenge in the Pure Deletion method,
one straightforward solution is to mask the deleted vertices as “deleted” and
push it into $Y$ in Algorithm 1. When we perform a greedy searching, the
“deleted” vertices are still visited. However, we do not push the “deleted”
vertices into the result priority queue—the visited “deleted” vertices are not
counted as top-$K$ candidates. The connections of the vertices in the graph
are unchanged. We illustrate this method in Figure 1 (c).

One advantage of the masking method is that the connectivity of the original
proximity graph is preserved and the expired vertices are still useful for
guiding the greedy search. However, the disadvantages are still obvious.
First, after masking, the search space for each query is large and
unnecessary. The visit to expired vertices slows down the greedy search speed.
Second, as all expired vertices are still stored, space grows continuously as
the stream performs, which may cause inevitable memory issues.

1: Input: deleted vertex $x_{i}$, proximity graph $G$ and its reverse graph
$G^{\prime}$, measurement function $f$, maximum outgoing degree of graph $d$.

2: $N^{\prime}(x_{i})\leftarrow$ out-neighbors of $x_{i}$ in $G^{\prime}$.

3: $N(x_{i})\leftarrow$ out-neighbors of $x_{i}$ in $G$.

4: for $x_{j}$ in $N(x_{i})$ do

5: Remove $(x_{i},x_{j})$ in $G$ and $(x_{j},x_{i})$ in $G^{\prime}$

6: end for

7: for $x_{k}$ in $N^{\prime}(x_{i})$ do

8: Remove $(x_{k},x_{i})$ in $G$ and $(x_{i},x_{k})$ in $G^{\prime}$

9: end for

10: Output: $G$, $G^{\prime}$.

Algorithm 4 PURE-DELETE$(x_{i},G,G^{\prime},f,d)$

1: Input: deleted vertex $x_{i}$, proximity graph $G$ and its reverse graph
$G^{\prime}$, measurement function $f$, maximum outgoing degree of graph $d$.

2: $N^{\prime}(x_{i})\leftarrow$ out-neighbors of $x_{i}$ in $G^{\prime}$.

3: $N(x_{i})\leftarrow$ out-neighbors of $x_{i}$ in $G$.

4: for $x_{j}$ in $N^{\prime}(x_{i})$ do

5: $N(x_{j})\leftarrow$ out-neighbors of $x_{j}$ in $G$.

6: $z\leftarrow\text{\small SELECT-
NEIGHBORS}(x_{j},N(x_{i}),1,N(x_{j})\cup\\{x_{j}\\})$.

7: Remove $(x_{j},x_{i})$ in $G$ and $(x_{i},x_{j})$ in $G^{\prime}$

8: if $z!=\mathsf{null}$ then

9: Add edges $(x_{j},z)$ to $G$ and $(z,x_{j})$ to $G^{\prime}$

10: end if

11: end for

12: Output: $G$, $G^{\prime}$.

Algorithm 5 LOCAL-RECONNECT$(x_{i},G,G^{\prime},f,d)$

###  5.3 Local Reconnect

To remedy the connectivity of the proximity graph while reducing the useless
visit in the search phase, another solution is to reset local connections
between the in-neighbors and out-neighbors of deleted vertex. The local
reconnect algorithm is presented in Algorithm 5 and Figure 1 (d). For each
expired vertex $x_{i}$, we first acquire its in-neighbors $N^{\prime}(x_{i})$
from $G^{\prime}$ and out-neighbors $N(x_{i})$ from $G$. Then, for each in-
neighbor $x_{j}$, we remove the edge from $x_{j}$ to $x_{i}$ in $G$ and
$G^{\prime}$. The last step is to select the most diverse vertex $z$ for
$x_{j}$ from $N(x_{i})$ and set out-going links from $x_{j}$ to $z$.

The idea of local reconnect is to compensate an edge for the in-neighbors of
deleted vertex. According to the nearest neighbor expansion, the out-neighbors
of expired vertex $x_{i}$ may be the neighbors of $x_{j}\in
N^{\prime}(x_{i})$. However, as the proximity graph is an approximation to
real Delaunay Graph, out-neighbors of $x_{i}$ may not be $x_{i}$’s exact
nearest neighbors. Therefore, the compensated edge may still not be good.

###  5.4 Global Reconnect

To compensate the in-neighbors of deleted vertices with better edges, we
propose a global reconnection algorithm for proximity graph update. The
algorithm is presented in Algorithm 6 and Figure 1 (e). For each expired
vertex $x_{i}$, we first acquire its in-neighbors $N^{\prime}(x_{i})$ from
$G^{\prime}$ and out-neighbors $N(x_{i})$ from $G$. After that, for every
element $x_{j}$ in $N^{\prime}(x_{i})$, we perform greedy search to obtain its
nearest neighbors $C$ on graph $G$. Then we use SELECT-NEIGHBORS function to
select the vertices set $N$ from $Z$ to fit the regulation for graph degree
$d$. Finally, we remove all the previous connections of $x_{j}$ and reset
edges by vertices in $N$.

The general idea behind this deletion method is to treat each in-neighbor
$x_{j}$ of the expired vertex as a new vertex and reset $x_{j}$’s connections
by re-inserting it on proximity graph. In this way, $x_{j}$ has an opportunity
to improve its connection by re-selecting vertices globally, Therefore, the
proximity graph may become a better approximation to Delaunay Graph. The
disadvantage is the updating time cost. however, it can be amortized by the
number of query operations. In real-world recommendation systems or search
engine, the number of queries(users) may be much larger than the number of
data(Ads, items) and one query may occur in high frequency. Thus, this
deletion method may outperform peers in online graph-based ANN search.

1: Input: deleted vertex $x_{i}$, proximity graph $G$ and its reverse graph
$G^{\prime}$, priority queue length $k$, measurement function $f$, out-
neighbor degree threshold $d$.

12: $N^{\prime}(x_{i})\leftarrow$ out-neighbors of $x_{i}$ in $G^{\prime}$.

3: for $x_{j}$ in $N^{\prime}(x_{i})$ do

4: $C\leftarrow\text{\small GREEDY-SEARCH}(x_{j},G,k,f)$.

5: $N\leftarrow\text{\small SELECT-NEIGHBORS}(x_{j},C,d,\\{x_{i}\\})$.

6: Remove $N(x_{j})$ in $G$ and $N^{\prime}(x_{j})$ in $G^{\prime}$

7: for $z\in N$ do

8: Add edges $(x_{j},z)$ to $G$ and $(z,x_{j})$ to $G^{\prime}$.

9: end for

10: end for

11: Output: $G$, $G^{\prime}$.

Algorithm 6 GLOBAL-RECONNECT$(x_{i},G,G^{\prime},k,f,d)$

##  6 Experiments

The objective of the experimental evaluation is to investigate the performance
of our update algorithms on real datasets. The performance includes both the
quality and efficiency of the retrieving top-$K$ candidates and the execution
time for indexing and searching. Specifically, we target to answer the
following questions:

  * •

Does the proposed global reconnect algorithm improve the performance of the
proximity graph? How does it compare with other graph update algorithms?

  * •

Does the proposed global reconnect algorithm robust to different update
patterns? How is the proposed algorithm and baselines perform in different
update patterns?

  * •

Does the proposed global reconnect algorithm reduce the total execution time
for insertion, deletion, and query?

Implementation. We implement the update algorithms as a C++11 prototype and
compile the code with g++-5.4.0 and “O3” optimization. The implementation in
Zhao et al. (2020) is employed as the proximity graph searching algorithm. The
code contains special functions to harness detailed profiling data.

Hardware System. We execute the experiments on a single node server. The
server has one Intel(R) Core(TM) i7-5960X CPU @ 3.00GHz (64 bit). It has 8
cores 16 threads, and 32 GB memory. Ubuntu 16.04.4 $64$-bit is the operating
system.

Retrieval Recall. In this paper, we use recall to measure the quality of
search, which is defined as the ratio of retrieved correct items over the
total correct items. The increase in recall relates to higher quality.

Data. We use $4$ real ANN benchmark datasets for experiments: SIFT (Jégou et
al., 2011b), GloVe200 (Pennington et al., 2014), NYTimes (Dua and Graff, 2017)
and GIST (Sandhawalia and Jégou, 2010). The dimension of vectors in these
datasets are 128,200,256 and 960. The NYTimes contains $280,000$ vectors while
others contain $1,000,000$ vectors. We have a diverge distribution of the
datasets—GloVe200 and NYTimes are more skewed compare to others.

Workload. For each of the 4 ANN datasets, we build 10 step workloads. Given
the base set, each step of the workload removes a set with 10,000 vectors and
then digests an set with 10,000 vectors. After that, a query set of 10,000
vectors are fed for top-K ANN search. In this case, the size for base set in
NYTimes is 180,000 while other sets are 900,000.

We consider two update patterns in the experiments: (a) Random updates: we
permute the whole data set before partition. Each vector of base, delete,
insert and query set is randomly selected from the original dataset. (b)
Cluster updates: we perform 10 class K-means clustering on the whole dataset.
Then we place the clusters in a sequence and build base, delete and insert
sets in the order of this sequence. Therefore, each step we remove several
clusters of vectors on graph. And then insert several clusters.

Methods. The 4 proposed algorithms are abbreviated as PURE, MASK, LOCAL and
GLOBAL. We also include ReBuild, which reconstructs the whole graph in each
update batch before query.

Figure 2: Relative queries per second to obtain 0.8 recall in each batch. The
update pattern for each batch is random updates.

###  6.1 Query Time

In this section, we present the query efficiency for our proposed algorithms
and baseline. In Figure 2 and Figure 3, we plot relative queries processed per
second (QPS) versus the query batches of workloads from 3 ANN datasets. The
queries processed per second (QPS) versus the query batches of workloads of
GIST dataset is also presented at Figure 3. In each workload, we follow the
pattern described above. Here the relative QPS is the result of real QPS
divided by the QPS of ReBuild algorithm in 0.8 Top-10, Top-20 or Top-100
recall.

Figure 3: Relative queries per second to obtain 0.8 recall in each batch. The
update pattern for each batch is clustered updates.

####  6.1.1 Random updates

The results in Figure 2 demonstrate that global reconnect achieves superior
performance than other delete-update-edge algorithms. For each dataset, the
QPSes of different methods start at the same point in query batch 0, which
means that they have the same performance when doing ANN search on the base
set without an update. Then, as the batch increases, we observe that our
global reconnect algorithm has higher QPS than other incremental update
methods. In SIFT dataset, the global reconnect algorithm even outperforms
reconstruct graph in QPS at 0.8 Top-10 recall. The same phenomenon also
happens in GloVe200 dataset and GIST dataset. For baseline methods, local
reconnect achieves the best performance. Pure deletion algorithm has a QPS
drop as batch increases, which means the proximity graph becomes inefficiency
due to the broken of connectivity. The vertex mask method achieves the lowest
QPS in every batch, illustrating that unnecessary visits of expired vertices
significantly influence the performance of ANN search on proximity graph.

####  6.1.2 Clustered updates

The results in Figure 3 demonstrates the superior performance of global
reconnect algorithm in workloads with clustered updates. Similar to Figure 2,
global reconnect achieves the best QPS in each batch at 0.8 Top-10, Top-20 or
Top-100 recall. In the clustered updates settings, we delete some vectors and
their nearest neighbors in each update step. Therefore, for vectors located
near to these clusters, a large fraction of their nearest neighbors are also
deleted. To update their connections, it requires global search on graph to
find new nearest neighbors.

Different from random updates, in clustered updates, we observe that our
global reconnect algorithm outperforms graph reconstruction in almost all
scenarios. Since the graph is constructed incrementally—new edges of an
inserted vertex are connected based on the existing vertices, those edges do
not accurately reveal the exact nearest neighbors. Our global reconnect
algorithm updates edges based on the entire graph. The edges in the graph are
refined during the updates—thus, a better performance than the reconstruction
is observed.

###  6.2 Overall Execution Time

To evaluate the execution time of our proposed global reconnect algorithm and
other baselines, we plot the accumulated time versus the number of operations.
The accumulated time grows as the number of operations increases. Lower plot
in the figure represents less workload execution time and better performance.
In this experiments, we consider two major characteristics of real-world
scenarios. First, embedding vectors of different usernames, are a much larger
group than the online product embeddings. Therefore, the number of queries may
be greater or equal to the number of vectors in database. Second, one single
user may look for products several times within an hour. Thus, one query
vector may be executed hundreds of times. Combining this two features, we
modify the original workload for SIFT. Given 900k base vectors, in every step
of the workload, we remove 10k vectors and insert the same amount of vectors,
but query with 200k, 1 million or 20 million vectors. The query set is
obtained by duplicating original 10k query set. Plots for these 3 workloads
corresponds to sub-figures from left to right in Figure 4.

Figure 4: Total execution time vs number of operations

In Figure 4 on the left, we show that if we query 200k vectors per batch.
Graph reconstruction has the highest time cost because the time cost
introduced by reset all vertices on graph. Then, the vertex masking is also
time inefficient. Although it deletes vertices in the fastest way, the
significant QPS drop in query phase makes it time-consuming. The global
reconnect, local reconnect and pure deletion have similar total execution
time. Global reconnect and local reconnect have higher execution time due to
the cost introduced by reconnection. Then, as the number of queries per batch
increases, we observe that our global reconnect algorithm achieves the lowest
total execution time.

This phenomenon validates the assumption that the updating time cost can be
amortized by the time saved by efficient ANN search. The same compensation
effect is also observed in graph reconstruction and local reconnect. As the
number of queries increases each batch, the total execution time of graph
reconstruction outperforms vertex masking and becomes closer to pure deletion.
The performance of local reconnect also improves as the margin between it and
pure deletion increases when we have more queries. Therefore, in real-world
applications that require online ANN search. Given the large amount of query
requests, the proposed IPGM framework outperforms static proximity graph based
methods due to the advantages taken in each search phase. Meanwhile, IPGM
reduces the latency of data deletion by avoiding reindexing the whole dataset
as a graph in each step.

###  6.3 Further Discussions

  * •

To obtain a level of recall at 0.8, the proposed global reconnect algorithm
improves the performance of incremental proximity graph in online ANN search
by a factor of 10.6$\%$ on average. In some scenarios, global reconnect
algorithm also outperforms graph reconstruction by a maximum factor of
18.8$\%$.

  * •

The proposed global reconnect algorithm is robust to random update pattern and
clustered update pattern. In clustered update when we delete a vector and its
nearest neighbors together, the global reconnect algorithm is a better way to
maintain a good approximation to Delaunay Graph.

  * •

Our proposed global reconnect algorithm reduces the total execution time when
the number of queries increases. In real-world settings where the number of
queries and data vectors are at the same level, our method has the fastest
execution time compared to other methods.

##  7 Conclusion

Graph-based ANN search methods index the dataset as proximity graph and has
shown superiority over other ANN search approaches. However, the inefficiency
of proximity graph against incremental data vector insertion and deletion
prevents it to be used in large-scale online ANN search. This paper presents a
proximity graph maintenance algorithm for online ANN search. We define the
incremental proximity graph update problem formally and give theoretical
proofs to demonstrate the feasibility of incremental vertex insertion and
deletion on the proximity graph. We propose 4 proximity graph update
algorithms to preserve the advantages of proximity graph in ANN search when
vertices are removed from the graph in different patterns. The experiments on
4 public ANN datasets show that our proposed algorithms outperform the
reconstruction baseline by a maximum factor of 18.8$\%$. The proposed updating
scheme eliminates the performance drop in online ANN methods on proximity
graphs, which enables the deployment in real-world search engines and
recommendation systems.

## References

  * André et al. (2015) Fabien André, Anne-Marie Kermarrec, and Nicolas Le Scouarnec.  Cache locality is not enough: High-performance nearest neighbor search with product quantization fast scan.  _Proc. VLDB Endow._ , 9(4):288–299, 2015.
  * Beaumont et al. (2007) Olivier Beaumont, Anne-Marie Kermarrec, and Etienne Riviere.  Peer to peer multidimensional overlays: Approximating complex structures.  In _Proceedings of Principles of Distributed Systems, 11th International Conference (OPODIS)_ , pages 315–328, Guadeloupe, French West Indies, 2007.
  * Berchtold et al. (1996) Stefan Berchtold, Daniel A. Keim, and Hans-Peter Kriegel.  The X-tree : An index structure for high-dimensional data.  In _Proceedings of 22th International Conference on Very Large Data Bases (VLDB)_ , pages 28–39, Mumbai (Bombay), India, 1996.
  * Broder (1997) Andrei Z. Broder.  On the resemblance and containment of documents.  In _Proceedings of the Conference on Compression and Complexity of SEQUENCES_ , pages 21–29, Positano, Amalfitan Coast, Salerno, Italy, 1997\.
  * Broder et al. (1998) Andrei Z. Broder, Moses Charikar, Alan M. Frieze, and Michael Mitzenmacher.  Min-wise independent permutations.  In _Proceedings of the Thirtieth Annual ACM Symposium on the Theory of Computing (STOC)_ , pages 327–336, Dallas, TX, 1998.
  * Cakir et al. (2017) Fatih Cakir, Sarah Adel Bargal, and Stan Sclaroff.  Online supervised hashing.  _Computer Vision and Image Understanding_ , 156:162–173, 2017.
  * Cao et al. (2016) Yue Cao, Mingsheng Long, Jianmin Wang, Han Zhu, and Qingfu Wen.  Deep quantization network for efficient image retrieval.  In _Proceedings of the Thirtieth AAAI Conference on Artificial Intelligence (AAAI)_ , pages 3457–3463, Phoenix, AZ, 2016.
  * Charikar (2002) Moses S. Charikar.  Similarity estimation techniques from rounding algorithms.  In _Proceedings on 34th Annual ACM Symposium on Theory of Computing (STOC)_ , pages 380–388, Montreal, Canada, 2002.
  * Datar et al. (2004) Mayur Datar, Nicole Immorlica, Piotr Indyk, and Vahab S. Mirrokni.  Locality-sensitive hashing scheme based on p-stable distributions.  In _Proceedings of the 20th ACM Symposium on Computational Geometry (SoCG)_ , pages 253–262, Brooklyn, NY, 2004.
  * Dave and Varma (2010) Kushal S. Dave and Vasudeva Varma.  Learning the click-through rate for rare/new ads from similar ads.  In _Proceeding of the 33rd International ACM SIGIR Conference on Research and Development in Information Retrieval (SIGIR)_ , pages 897–898, Geneva, Switzerland, 2010.
  * Dua and Graff (2017) Dheeru Dua and Casey Graff.  UCI machine learning repository, 2017.  URL http://archive.ics.uci.edu/ml.
  * Fan et al. (2019) Miao Fan, Jiacheng Guo, Shuai Zhu, Shuo Miao, Mingming Sun, and Ping Li.  MOBIUS: towards the next generation of query-ad matching in baidu’s sponsored search.  In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining (KDD)_, pages 2509–2517, Anchorage, AK, 2019.
  * Friedman et al. (1975) Jerome H. Friedman, F. Baskett, and L. Shustek.  An algorithm for finding nearest neighbors.  _IEEE Transactions on Computers_ , 24:1000–1006, 1975.
  * Friedman et al. (1977) Jerome H. Friedman, J. Bentley, and R. Finkel.  An algorithm for finding best matches in logarithmic expected time.  _ACM Transactions on Mathematical Software_ , 3:209–226, 1977.
  * Fu et al. (2019) Cong Fu, Chao Xiang, Changxu Wang, and Deng Cai.  Fast approximate nearest neighbor search with the navigating spreading-out graph.  _Proc. VLDB Endow._ , 12(5):461–474, 2019.
  * Gao et al. (2015) Jinyang Gao, H. V. Jagadish, Beng Chin Ooi, and Sheng Wang.  Selective hashing: Closing the gap between radius search and k-nn search.  In _Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD)_ , pages 349–358, Sydney, Australia, 2015.
  * Gao et al. (2020) Weihao Gao, Xiangjun Fan, Chong Wang, Jiankai Sun, Kai Jia, Wenzhi Xiao, Ruofan Ding, Xingyan Bin, Hui Yang, and Xiaobing Liu.  Deep retrieval: Learning a retrievable structure for large-scale recommendations.  _arXiv preprint arXiv:2007.07203_ , 2020.
  * Ge et al. (2013) Tiezheng Ge, Kaiming He, Qifa Ke, and Jian Sun.  Optimized product quantization for approximate nearest neighbor search.  In _Proceedings of the 2013 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 2946–2953, Portland, OR, 2013.
  * Ghashami and Abdullah (2015) Mina Ghashami and Amirali Abdullah.  Binary coding in stream.  _arXiv preprint arXiv:1503.06271_ , 2015.
  * Gionis et al. (1999) Aristides Gionis, Piotr Indyk, and Rajeev Motwani.  Similarity search in high dimensions via hashing.  In _Proceedings of 25th International Conference on Very Large Data Bases (VLDB)_ , pages 518–529, Edinburgh, Scotland, UK, 1999.
  * Guare (1990) John Guare.  _Six degrees of separation: A play_.  Vintage, 1990.
  * Hajebi et al. (2011) Kiana Hajebi, Yasin Abbasi-Yadkori, Hossein Shahbazi, and Hong Zhang.  Fast approximate nearest-neighbor search with k-nearest neighbor graph.  In _Proceedings of the 22nd International Joint Conference on Artificial Intelligence (IJCAI)_ , pages 1312–1317, Barcelona, Catalonia, Spain, 2011.
  * Indyk and Motwani (1998) Piotr Indyk and Rajeev Motwani.  Approximate nearest neighbors: Towards removing the curse of dimensionality.  In _Proceedings of the Thirtieth Annual ACM Symposium on the Theory of Computing (STOC)_ , pages 604–613, Dallas, TX, 1998.
  * Jégou et al. (2011a) Hervé Jégou, Matthijs Douze, and Cordelia Schmid.  Product quantization for nearest neighbor search.  _IEEE Trans. Pattern Anal. Mach. Intell._ , 33(1):117–128, 2011a.
  * Jégou et al. (2011b) Hervé Jégou, Romain Tavenard, Matthijs Douze, and Laurent Amsaleg.  Searching in one billion vectors: Re-rank with source coding.  In _Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP)_ , pages 861–864, Prague, Czech Republic, 2011b.
  * Karinthy (1929) Frigyes Karinthy.  Chain-links.  _Everything is different_ , pages 21–26, 1929.
  * Kleinberg (2000) Jon M. Kleinberg.  The small-world phenomenon: an algorithmic perspective.  In _Proceedings of the Thirty-Second Annual ACM Symposium on Theory of Computing (STOC)_ , pages 163–170, Portland, OR, 2000.
  * Leng et al. (2015) Cong Leng, Jiaxiang Wu, Jian Cheng, Xiao Bai, and Hanqing Lu.  Online sketching hashing.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 2503–2511, Boston, MA, 2015.
  * Li (2018) Ping Li.  Several tunable GMM kernels.  _arXiv preprint arXiv:1805.02830_ , 2018.
  * Li (2019) Ping Li.  Sign-full random projections.  In _Proceedings of the Thirty-Third AAAI Conference on Artificial Intelligence (AAAI)_ , pages 4205–4212, Honolulu, HI, 2019.
  * Li and Church (2005) Ping Li and Kenneth Ward Church.  Using sketches to estimate associations.  In _Proceedings of the Human Language Technology Conference and the Conference on Empirical Methods in Natural Language Processing (HLT/EMNLP)_ , pages 708–715, Vancouver, Canada, 2005.
  * Li and Zhang (2017) Ping Li and Cun-Hui Zhang.  Theory of the GMM kernel.  In _Proceedings of the 26th International Conference on World Wide Web (WWW)_ , pages 1053–1062, Perth, Australia, 2017.
  * Li et al. (2021) Ping Li, Xiaoyun Li, Gennady Samorodnitsky, and Weijie Zhao.  Consistent sampling through extremal process.  In _Proceedings of the Web Conference (WWW)_ , pages 1317–1327, 2021\.
  * Li and Li (2019) Xiaoyun Li and Ping Li.  Random projections with asymmetric quantization.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 10857–10866, Vancouver, Canada, 2019.
  * Lin et al. (1994) King-Ip Lin, H. V. Jagadish, and Christos Faloutsos.  The TV-tree: An index structure for high-dimensional data.  _VLDB J._ , 3(4):517–542, 1994.
  * Liong et al. (2015) Venice Erin Liong, Jiwen Lu, Gang Wang, Pierre Moulin, and Jie Zhou.  Deep hashing for compact binary codes learning.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 2475–2483, Boston, MA, 2015.
  * Malkov and Yashunin (2020) Yury Malkov and Dmitry Yashunin.  Efficient and robust approximate nearest neighbor search using hierarchical navigable small world graphs.  _IEEE Trans. Pattern Anal. Mach. Intell._ , 42(4):824–836, 2020.
  * Malkov et al. (2014) Yury Malkov, Alexander Ponomarenko, Andrey Logvinov, and Vladimir Krylov.  Approximate nearest neighbor algorithm based on navigable small world graphs.  _Inf. Syst._ , 45:61–68, 2014.
  * Morozov and Babenko (2018) Stanislav Morozov and Artem Babenko.  Non-metric similarity graphs for maximum inner product search.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 4726–4735, Montreal, Canada, 2018.
  * Navarro (2002) Gonzalo Navarro.  Searching in metric spaces by spatial approximation.  _The VLDB Journal_ , 11(1):28–46, 2002.
  * Paredes and Chávez (2005) Rodrigo Paredes and Edgar Chávez.  Using the _k_ -nearest neighbor graph for proximity searching in metric spaces.  In _Proceedings of the 12th International Conference on String Processing and Information Retrieval (SPIRE)_ , pages 127–138, Buenos Aires, Argentina, 2005.
  * Pazzani and Billsus (2007) Michael J. Pazzani and Daniel Billsus.  Content-based recommendation systems.  In _Proceedings of the Conference on the Adaptive Web, Methods and Strategies of Web Personalization_ , pages 325–341, 2007.
  * Pennington et al. (2014) Jeffrey Pennington, Richard Socher, and Christopher D. Manning.  Glove: Global vectors for word representation.  In _Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 1532–1543, Doha, Qatar, 2014.
  * Ram and Sinha (2019) Parikshit Ram and Kaushik Sinha.  Revisiting kd-tree for nearest neighbor search.  In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining (KDD)_, pages 1378–1388, Anchorage, AK, 2019.
  * Sakurai et al. (2000) Yasushi Sakurai, Masatoshi Yoshikawa, Shunsuke Uemura, and Haruhiko Kojima.  The a-tree: An index structure for high-dimensional spaces using relative approximation.  In _Proceedings of 26th International Conference on Very Large Data Bases (VLDB)_ , pages 516–526, Cairo, Egypt, 2000.
  * Sandhawalia and Jégou (2010) Harsimrat Sandhawalia and Hervé Jégou.  Searching with expectations.  In _Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP)_ , pages 1242–1245, Dallas, TX, 2010.
  * Shen et al. (2015) Fumin Shen, Chunhua Shen, Wei Liu, and Heng Tao Shen.  Supervised discrete hashing.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 37–45, Boston, MA, 2015.
  * Shrivastava and Li (2012) Anshumali Shrivastava and Ping Li.  Fast near neighbor search in high-dimensional binary data.  In _Proceedings of the European Conference on Machine Learning and Knowledge Discovery in Databases (ECML-PKDD), Part I_ , pages 474–489, Bristol, UK, 2012.
  * Singh et al. (2021) Aditi Singh, Suhas Jayaram Subramanya, Ravishankar Krishnaswamy, and Harsha Vardhan Simhadri.  FreshDiskANN: A fast and accurate graph-based ann index for streaming similarity search.  _arXiv preprint arXiv:2105.09613_ , 2021.
  * Tan et al. (2019) Shulong Tan, Zhixin Zhou, Zhaozhuo Xu, and Ping Li.  On efficient retrieval of top similarity vectors.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ , pages 5235–5245, Hong Kong, China, 2019.
  * Tan et al. (2020a) Shulong Tan, Zhixin Zhou, Zhaozhuo Xu, and Ping Li.  Fast item ranking under neural network based measures.  In _Proceedings of the Thirteenth ACM International Conference on Web Search and Data Mining (WSDM)_ , Houston, TX, 2020a.
  * Tan et al. (2020b) Shulong Tan, Zhixin Zhou, Zhaozhuo Xu, and Ping Li.  Fast item ranking under neural network based measures.  In _Proceedings of the Thirteenth ACM International Conference on Web Search and Data Mining (WSDM)_ , pages 591–599, Houston, TX, 2020b.
  * Tan et al. (2021a) Shulong Tan, Zhaozhuo Xu, Weijie Zhao, Hongliang Fei, Zhixin Zhou, and Ping Li.  Norm adjusted proximity graph for fast inner product retrieval.  In _Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery and Data Mining (KDD)_ , pages 1552–1560, Virtual Event, Singapore, 2021a.
  * Tan et al. (2021b) Shulong Tan, Weijie Zhao, and Ping Li.  Fast neural ranking on bipartite graph indices.  _Proc. VLDB Endow._ , 15(4):794–803, 2021b.
  * Vanderkam et al. (2013) Dan Vanderkam, Rob Schonberger, Henry Rowley, and Sanjiv Kumar.  Nearest neighbor search in google correlate.  Technical report, Google, 2013.
  * Wei et al. (2014) Ying Wei, Yangqiu Song, Yi Zhen, Bo Liu, and Qiang Yang.  Scalable heterogeneous translated hashing.  In _Proceedings of the 20th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD)_ , pages 791–800, New York, NY, 2014.
  * Weiss et al. (2008) Yair Weiss, Antonio Torralba, and Robert Fergus.  Spectral hashing.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 1753–1760, Vancouver, Canada, 2008.
  * Xu et al. (2018a) Donna Xu, Ivor W. Tsang, and Ying Zhang.  Online product quantization.  _IEEE Trans. Knowl. Data Eng._ , 30(11):2185–2198, 2018a.
  * Xu et al. (2018b) Jun Xu, Xiangnan He, and Hang Li.  Deep learning for matching in search and recommendation.  In _Proceedings of the 41st International ACM SIGIR Conference on Research & Development in Information Retrieval (SIGIR)_, pages 1365–1368, Ann Arbor, MI, 2018b.
  * Yu et al. (2022) Tan Yu, Jie Liu, Yi Yang, Yi Li, Hongliang Fei, and Ping Li.  EGM: Enhanced graph-based model for large-scale video advertisement search.  In _Proceedings of the 28th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining (KDD)_, Washington, DC, 2022\.
  * Zezula et al. (1998) Pavel Zezula, Pasquale Savino, Giuseppe Amato, and Fausto Rabitti.  Approximate similarity retrieval with M-trees.  _VLDB J._ , 7(4):275–293, 1998.
  * Zhang et al. (2014) Ting Zhang, Chao Du, and Jingdong Wang.  Composite quantization for approximate nearest neighbor search.  In _Proceedings of the 31th International Conference on Machine Learning (ICML)_ , pages 838–846, Beijing, China, 2014.
  * Zhao et al. (2020) Weijie Zhao, Shulong Tan, and Ping Li.  SONG: approximate nearest neighbor search on GPU.  In _Proceedings of the 36th IEEE International Conference on Data Engineering (ICDE)_ , pages 1033–1044, Dallas, TX, 2020.
  * Zhou et al. (2019) Zhixin Zhou, Shulong Tan, Zhaozhuo Xu, and Ping Li.  Möbius transformation for fast inner product search on graph.  In _Proceedings of Advances in Neural Information Processing Systems (NeurIPS)_ , pages 8216–8227, Vancouver, Canada, 2019.
  * Zhu et al. (2016) Han Zhu, Mingsheng Long, Jianmin Wang, and Yue Cao.  Deep hashing network for efficient similarity retrieval.  In _Proceedings of the Thirtieth AAAI Conference on Artificial Intelligence (AAAI)_ , pages 2415–2421, Phoenix, AZ, 2016.
  * Zhu et al. (2018) Han Zhu, Xiang Li, Pengye Zhang, Guozheng Li, Jie He, Han Li, and Kun Gai.  Learning tree-based deep model for recommender systems.  In _Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining (KDD)_, pages 1079–1088, London, UK, 2018.
  * Zhu et al. (2019) Han Zhu, Daqing Chang, Ziru Xu, Pengye Zhang, Xiang Li, Jie He, Han Li, Jian Xu, and Kun Gai.  Joint optimization of tree-based index and deep model for recommender systems.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 3973–3982, Vancouver, Canada, 2019.
  * Zhuo et al. (2020) Jingwei Zhuo, Ziru Xu, Wei Dai, Han Zhu, Han Li, Jian Xu, and Kun Gai.  Learning optimal tree models under beam search.  In _Proceedings of the 37th International Conference on Machine Learning (ICML)_ , pages 11650–11659, Virtual Event, 2020.


# Statistical anisotropy of primordial gravitational waves from generalized
$\delta N$ formalism

Takahiro Tanaka<EMAIL_ADDRESS>Department of Physics,
Kyoto University, Kyoto 606-8502, Japan  Center for Gravitational Physics and
Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto
University, Kyoto 606-8502, Japan  Yuko Urakawa<EMAIL_ADDRESS>Institute
of Particle and Nuclear Studies, High Energy Accelerator Research Organization
(KEK), Oho 1-1, Tsukuba 305-0801, Japan  The Graduate University for Advanced
Studies (SOKENDAI), Tsukuba 305-0801, Japan

###### Abstract

In this letter, we demonstrate how to use the generalized $\delta N$
formalism, which enables us to compute all the large scale fluctuations,
including the gravitational waves, solely by solving the evolution of the
background homogeneous Universe. Using the Noether charge density, we derive
an analytic formula which describes the mapping between the fluctuations at
the horizon crossing and the sourced gravitational waves at the end of
inflation. This formula can apply also to an inflation model with an
anisotropic background.

Introduction.– The $\delta N$ formalism [1, 2, 3], which is based on the
separate universe approach [4, 5], has played the central role to connect the
prediction of each inflation model with various observations, enabling a
simple computation of the superhorizon dynamics and also providing an
intuitive understanding on the evolution of the primordial fluctuations.
However, their application was limited to a system which contains only scalar
fields and it could not be used even to compute the gravitational waves. In
Ref. [6], we have shown that the separate universe approach can be applied
under a rather general setup. Based on this, we can generalize the $\delta N$
formalism to compute the adiabatic perturbation and the gravitational waves
sourced by non-zero spin fields. In this letter, we demonstrate this
computation, considering a model with $U(1)$ gauge fields.

g$\delta N$ formalism.– The $\delta N$ formalism computes the evolution of
fluctuations at the leading order of the gradient expansion [4, 7], which is
an expansion scheme with respect to the spatial gradient. The gradient
expansion starts with smoothing the small scale fluctuations. As a consequence
of the smoothing, operating the spatial gradient gives rise to the suppression
by a small quantity $\epsilon$, which is usually characterized by the spatial
variation of $\varphi^{a}$ within each causally connected patch. At the
leading order of the gradient expansion, we simply send $\epsilon$ to 0. As
shown in [6], the separate universe approach and the $\delta N$ formalism
generically applies to a theory which satisfies the spatial diffeomorphism
invariance and the locality.

Under these two conditions, one can compute the time evolution of the
inhomogeneous Universe simply by solving a set of the corresponding ordinary
differential equations for different initial conditions specified at around
the horizon crossing. The conventional $\delta N$ formalism can be applied
only to a system with scalar fields, while the generalized $\delta N$
formalism (g$\delta N$ formalism) can be applied broadly to a general model
that satisfies the above-mentioned two conditions [6].

Preliminaries.– In this letter, we consider a system with $D$ scalar fields
$\phi^{I}$ and $D^{\prime}$ U(1) gauge fields $A^{\mu}_{(\alpha)}$ whose
Lagrangian density is given by

$\displaystyle{\cal L}_{\rm mat}$ $\displaystyle={\cal
P}(X^{IJ},\,\phi^{I})-\frac{f_{(\alpha)}^{2}(\phi^{I})}{4}F_{\mu\nu(\alpha)}F^{\mu\nu}_{(\alpha)}\,,$
(1)

where $F_{\mu\nu(\alpha)}$ denotes the field strengths of the gauge fields
$A^{\mu}_{(\alpha)}$. Here and hereafter, we do not explicitly write the
summation over $\alpha$, which counts different gauge fields.

In this letter, we express the 4-dimensional line element as

$\displaystyle ds^{2}=-N^{2}dt^{2}+g_{ij}(dx^{i}+N^{i}dt)(dx^{j}+N^{j}dt)\,,$
(2)

with $i,j=1,\cdots,3$. We express the spatial metric as
$g_{ij}=e^{2\psi}\,\gamma_{ij}$, where $\gamma_{ij}$ satisfies
$\det[\gamma]=1$. Using $\psi$, the determinant of $g_{ij}$ is given by
$g=e^{6\psi}$. In the separate universe evolution, we need to employ local
gauge conditions [6]. Here, we adopt $N_{i}=0$ and $A_{0}=0$. In $N_{i}=0$
gauge, the expansion $K$ and the shear ${A^{i}}_{j}$, given by the trace part
and traceless part of the extrinsic curvature, read $K=3\dot{\psi}/N$ and
${A^{i}}_{j}=\gamma^{im}\dot{\gamma}_{mj}/2N$, respectively, with a dot being
the time derivative with respect to $t$. We determine the time slicing at the
horizon crossing, requiring $\delta\psi(t_{*},\,\hbox{\boldmath{$x$}})=0$ and
the one at $t=t_{\rm f}$, requiring $\delta K(t_{\rm
f},\,\hbox{\boldmath{$x$}})=0$. The residual gauge degrees of freedom are
eliminated by imposing additional gauge conditions on the initial condition at
the horizon crossing [8]. We introduce $\delta{\gamma^{i}}_{j}$ and
$\delta\gamma_{ij}(t,\,\hbox{\boldmath{$x$}})$ as
$\gamma_{ij}(t,\,\hbox{\boldmath{$x$}})\equiv\bar{\gamma}_{ik}(t)\left[e^{\delta\gamma(t,\,\hbox{\boldmath{\scriptsize$x$}})}\right]^{k}\\!_{j}$
and
$\delta\gamma_{ij}(t,\,\hbox{\boldmath{$x$}})\equiv\bar{\gamma}_{ik}(t)\delta{\gamma^{k}}_{j}(t,\,\hbox{\boldmath{$x$}})$.
Here and hereafter, we put a bar on background variables.

The Maxwell equations are given by $\partial_{t}\pi_{(\alpha)}^{i}={\cal
O}(\epsilon^{2})$ with the conjugate momenta being

$\displaystyle\pi_{(\alpha)}^{i}\equiv\frac{\partial(N\sqrt{g}{\cal L}_{\rm
matter})}{\partial\dot{A}_{i(\alpha)}}=\sqrt{g}f_{(\alpha)}^{2}(\phi^{I})g^{ij}\frac{\dot{A}_{j(\alpha)}}{N}\,,$

which correspond to the Noether charge densities at the leading order of the
gradient expansion [6]. The energy densities of the gauge fields are given by

$\displaystyle\rho_{A(\alpha)}=\frac{\gamma_{ij}\pi_{(\alpha)}^{i}\pi_{(\alpha)}^{j}}{2f_{(\alpha)}^{2}e^{4\psi}}+{\cal
O}(\epsilon^{2})\,,$ (3)

which indicates that $\rho_{A(\alpha)}$ grows, when $f_{(\alpha)}$ decreases
faster than $e^{-2\psi}$.

Primordial gravitational waves.– In this letter, we assume the near de Sitter
evolution (but not the slow-roll approximation) and that the shear is much
smaller than the expansion, i.e.,

$\displaystyle({\rm A1})\quad-\frac{3\dot{K}}{NK^{2}}\ll
1\,,\quad\frac{{A^{i}}_{j}{A^{j}}_{i}}{K^{2}}\ll 1\,.$ (4)

Using the solution of the Maxwell equations and discarding trivial decaying
solutions, we can formally solve the traceless part of $(i,\,j)$ component of
the Einstein equations as

$\displaystyle\gamma_{ij}(t,\,\hbox{\boldmath{$x$}})=\gamma_{ij*}(\hbox{\boldmath{$x$}})-2\left[\gamma_{il*}(\hbox{\boldmath{$x$}})\gamma_{jm*}(\hbox{\boldmath{$x$}})\right]^{\rm
TL}$
$\displaystyle\qquad\qquad\quad\times\hat{\pi}_{(\alpha)}^{l}(\hbox{\boldmath{$x$}})\hat{\pi}_{(\alpha)}^{m}(\hbox{\boldmath{$x$}})\,\hat{{\cal
I}}_{(\alpha)}(t;\\{\varphi^{a^{\prime}}_{*}\\}^{\prime})+{\cal
O}(\epsilon^{2})\,,$ (5)

with $\hat{{\cal I}}_{(\alpha)}$ and $\hat{\pi}^{i}_{(\alpha)}$ being

$\displaystyle\hat{{\cal
I}}_{(\alpha)}(t;\\{\varphi^{a^{\prime}}_{*}\\}^{\prime})$
$\displaystyle\equiv\frac{\gamma_{ij*}\pi_{(\alpha)}^{i}\pi^{j}_{(\alpha)}}{M_{\rm
pl}^{2}}\\!\int^{t}_{t_{*}}\\!\frac{dt^{\prime}N}{e^{3\psi(t^{\prime})}}\\!\int^{t^{\prime}}_{t_{*}}\\!\frac{dt^{\prime\prime}N}{e^{\psi(t^{\prime\prime})}f_{(\alpha)}^{2}(\phi^{I}(t^{\prime\prime}))}\,,$

and

$\displaystyle\hat{\pi}_{(\alpha)}^{i}\equiv\frac{\pi_{(\alpha)}^{i}}{\sqrt{\gamma_{ij*}\pi_{(\alpha)}^{i}\pi_{(\alpha)}^{j}}}=\frac{\gamma^{ij}_{*}\dot{A}_{j(\alpha)*}}{\sqrt{\gamma^{kl}_{*}\dot{A}_{k(\alpha)*}\dot{A}_{l(\alpha)*}}}\,,$
(6)

where $\\{\varphi^{a^{\prime}}_{*}\\}^{\prime}$ denotes the set of the fields
by which we provide the initial condition at the horizon crossing (a detailed
argument can be found in [8]). Here, TL denotes the traceless part about the
$(i,j)$ indices defined by using the spatial metric $\gamma_{ij*}$. The formal
solution (5) can be more easily obtained by using the Noether charge density
[6, 8] than by solving the Einstein equations directly. The second assumption
in (A1) can be verified when the energy density of the gauge fields,
$\rho_{A}=\sum_{\alpha}\rho_{A(\alpha)}$, remains much smaller than the total
one. Here we only take into account the leading contribution of the gauge
fields to $\gamma_{ij}$. This equation also gives the longitudinal part of
$\gamma_{ij}$, which is necessary to compute $\zeta$ in the present gauge.
Since $\hat{\pi}_{(\alpha)}^{i}$ is constant regardless the evolution of
$f_{(\alpha)}$, the scalar fields contribute to $\gamma_{ij}$ only through
$\hat{{\cal I}}_{(\alpha)}$.

Ignoring the time variation of the expansion under the assumption (4), we can
further rewrite $\hat{\cal I}_{(\alpha)f}\equiv\hat{\cal I}_{(\alpha)}(t_{f})$
with $t_{f}$ being the reheating surface as

$\displaystyle\hat{\cal
I}_{(\alpha)f}\simeq\frac{1}{\rho_{*}}\int^{\psi_{f}}_{\psi_{*}}\frac{d\psi^{\prime\prime}\gamma_{ij}\pi_{(\alpha)}^{i}\pi_{(\alpha)}^{j}}{e^{4\psi^{\prime\prime}}f_{(\alpha)}^{2}(\phi^{I}(\psi^{\prime\prime}))}$
(7)

with $\rho_{*}$ being the total energy density at $t_{*}$. Here, we inserted
$\gamma_{ij*}\pi^{i}\pi^{j}$ in the integral, since the time variation of
$\gamma_{ij}$ becomes higher order under the assumption (A1). Here, we have
integrated over $t^{\prime}$ after exchanging the order of the integrals and
ignored the exponentially suppressed term. Using Eq. (3), one can find that
the integrand of $\hat{\cal I}_{(\alpha)f}$ corresponds to the energy density
of each gauge field. Therefore, when $\rho_{A(\alpha)}$ takes a maximum value
$\rho_{A(\alpha)}^{\rm max}$ during $\Delta\psi_{(\alpha)}$ after the horizon
crossing, $\hat{\cal I}_{(\alpha)f}$ is roughly given by $\hat{\cal
I}_{(\alpha)f}\sim 2\Delta\psi_{(\alpha)}\rho_{A(\alpha)}^{\rm max}/\rho_{*}$.

We further assume

$\displaystyle({\rm A2})\quad\left|\frac{\delta\hat{\cal
I}_{(\alpha)f}}{\hat{\cal
I}_{(\alpha)f}}\right|\ll|\delta\hat{\pi}_{(\alpha)}^{i}|\,.$ (8)

For example, when $\rho_{A(\alpha)}^{\rm max}$ is determined by the potential
gradient $V_{I}/V$ as a consequence of the balance with the friction, the
fluctuation of $\rho_{A(\alpha)}^{\rm max}$ is suppressed by the slow-roll
parameter [9, 10]. Meanwhile, in the model discussed in Ref. [11], since
$\rho_{A(\alpha)}^{\rm max}$ becomes constant, there is no fluctuation. Then,
the left hand side of (A2) becomes
$|\delta\Delta\psi_{(\alpha)}|/\Delta\psi_{(\alpha)}$, which typically becomes
smaller than
$|\delta\hat{\pi}^{i}_{(\alpha)}|\\!\sim\\!|\delta\rho_{A(\alpha)*}/\rho_{A(\alpha)*}|$
by a factor of $1/\Delta\psi_{(\alpha)}$.

Perturbative expansion in g$\delta N$ formalism.– In the conventional $\delta
N$ formalism, the adiabatic curvature perturbation $\zeta$ is given by
computing the $e$-folding number under different initial conditions.
Similarly, we can compute the gravitational waves simply by computing the
anisotropic expansion. Under the assumptions listed above, perturbing Eq. (5),
the dominant contribution at the linear perturbation is given by

$\displaystyle\delta\gamma_{ij}(t_{f},\,\hbox{\boldmath{$x$}})$ (9)
$\displaystyle\simeq\delta\gamma_{ij*}(\hbox{\boldmath{$x$}})-4\left[\bar{\gamma}_{il*}\bar{\gamma}_{jm*}\right]^{\rm
TL}\bar{\hat{\pi}}^{l}_{(\alpha)}\delta\hat{\pi}^{m}_{(\alpha)}(\hbox{\boldmath{$x$}})\,\bar{\hat{{\cal
I}}}_{f(\alpha)}\,.$ (10)

The first term in the second line is the usual vacuum contribution. The second
term describes the shear fluctuation sourced by the fluctuations in the
directions of $\pi^{i}_{(\alpha)}$. The fluctuation of $\gamma_{ij*}$ also
perturbs the second term of Eq. (5), yielding several terms whose amplitudes
amount to $\bar{\hat{{\cal I}}}_{f(\alpha)}(H_{*}/M_{\rm
pl})\sim(\bar{\rho}^{\rm
max}_{A(\alpha)}/\bar{\rho}_{*})\Delta\psi_{(\alpha)}(H_{*}/M_{\rm pl})$. They
turn out to be smaller than the second term of Eq. (10) by
$\sqrt{\bar{\rho}^{\rm max}_{A(\alpha)}/\bar{\rho}_{*}}<1$. This suppression
originates from the difference between the amplitudes of $\delta\gamma_{ij*}$
and $\delta\hat{\pi}^{i}_{(\alpha)*}$.

Perturbing $\hat{\pi}^{i}_{(\alpha)}$, we obtain

$\displaystyle\delta\hat{\pi}_{(\alpha)}^{i}=\frac{\bar{\gamma}^{ij}_{*}\delta\dot{A}_{j(\alpha)*}}{\dot{\bar{A}}_{(\alpha)*}}-\frac{\bar{\gamma}_{*}^{kl}\dot{\bar{A}}_{k(\alpha)*}\delta\dot{A}_{l(\alpha)*}\bar{\gamma}^{ij}_{*}\dot{\bar{A}}_{j(\alpha)*}}{\dot{\bar{A}}_{(\alpha)*}^{3}}+\cdots\,,$

where we have introduced
$\dot{\bar{A}}_{(\alpha)*}\equiv\sqrt{\bar{\gamma}^{kl}_{*}\dot{\bar{A}}_{k(\alpha)*}\dot{\bar{A}}_{l(\alpha)*}}$.
Here we abbreviated the terms with the fluctuation of the spatial metric,
which only give sub-dominant contributions. Without these abbreviated terms,
we find
$\bar{\gamma}_{ij*}\delta\hat{\pi}_{(\alpha)}^{i}\bar{\pi}^{j}_{(\alpha)}=0$.

Polarization bases.–In what follows, we set the background spatial metric at
the reheating surface $t_{f}$ as $\bar{\gamma}_{ij}(t_{f})=\delta_{ij}$. Then,
at $t=t_{f}$, we can define the gravitational waves as usual by using the
polarization tensors $e_{ij}^{(\lambda_{\rm gw})}$ with $\lambda_{\rm
gw}=+,\,\times$, which satisfy $\hat{k}_{i}e_{ij}^{(\lambda_{\rm
gw})}(\hat{\hbox{\boldmath{$k$}}})=0$, $e_{ii}^{(\lambda_{\rm
gw})}(\hat{\hbox{\boldmath{$k$}}})=0$, and $e_{ij}^{(\lambda_{\rm
gw})}(\hat{\hbox{\boldmath{$k$}}})e_{ij}^{(\lambda^{\prime}_{\rm
gw})}(\hat{\hbox{\boldmath{$k$}}})=\delta^{\lambda_{\rm
gw}\lambda^{\prime}_{\rm gw}}$. The adiabatic curvature perturbation is also
given by the usual definition as
$\zeta(t_{f})=\delta\psi(t_{f})-\frac{1}{4}\hat{k}_{i}\hat{k}_{j}\delta\gamma_{ij}(t_{f})$.
Here, $\hat{\hbox{\boldmath{$k$}}}$ denotes the unit wavenumber
$\hbox{\boldmath{$k$}}/k$. Here and hereafter, we lower and raise the indices
$\hat{k}^{i}$ and $e_{ij}^{(\lambda_{\rm gw})}$ by using
$\bar{\gamma}_{ij}(t_{f})=\delta_{ij}$ and
$\bar{\gamma}^{ij}(t_{f})=\delta^{ij}$. (When the background shear had already
become negligible at $t=t_{f}$, the background spatial metric remains
$\bar{\gamma}_{ij}(t)=\delta_{ij}$ all the time after $t_{f}$, ensuring the
linear decomposition among the scalar, vector, and tensor type perturbations.)
Using the two polarization bases of the gauge fields,
$e^{(A)}_{i}(\hat{\hbox{\boldmath{$k$}}})$ with $A=1,\,2$, we define
$e_{ij}^{(\lambda_{\rm gw})}(\hat{\hbox{\boldmath{$k$}}})$ as
$e^{(+)}_{ij}\equiv(e_{i}^{(1)}e_{j}^{(1)}-e_{i}^{(2)}e_{j}^{(2)})/\sqrt{2}$
and
$e^{(\times)}_{ij}\equiv(e_{i}^{(1)}e_{j}^{(2)}+e_{i}^{(2)}e_{j}^{(1)})/\sqrt{2}$.
For an arbitrary $\hat{k}^{i}$, we can choose $e_{i}^{(2)}$ such that being
orthogonal to one of the background gauge fields, i.e.,
$e_{i}^{(2)}\bar{\pi}^{i}_{(1)}=0$. In fact, we can choose $e_{i}^{(1)}$ in
the 2D plane spanned by $\hat{k}_{i}$ and $\bar{\pi}_{(1)}^{i}$, ensuring
$e_{i}^{(2)}\bar{\pi}^{i}_{(1)}=0$. When we define $\Theta_{\alpha}$ as
$\cos\Theta_{\alpha}\equiv\hat{\hbox{\boldmath{$k$}}}\cdot\bar{\hat{\hbox{\boldmath{$\pi$}}}}_{(\alpha)}$,
under this choice, we obtain
$\hbox{\boldmath{$e$}}^{(1)}\cdot\bar{\hat{\hbox{\boldmath{$\pi$}}}}_{(1)}=\cos(\pi/2+\Theta_{1})=-\sin\Theta_{1}$.

Linear polarization.–For $D^{\prime}=1$, choosing the direction of
$e_{i}^{(A)}$ to satisfy $e_{i}^{(2)}\bar{\pi}^{i}=0$ and operating
$e_{ij}^{(\lambda_{\rm gw})}$ on Eq. (10), we obtain

$\displaystyle\gamma^{(\lambda_{\rm gw})}(t_{f},\,\hbox{\boldmath{$k$}})$
$\displaystyle\\!=\gamma^{(\lambda_{\rm
gw})}_{*}(\hbox{\boldmath{$k$}})+2\sqrt{2}\sin\Theta\,\bar{\hat{{\cal
I}}}_{f}\\!\left\\{\begin{array}[]{ll}\\!\\!\cos^{2}\Theta\frac{e^{(1)}_{i}\delta\dot{A}_{i*}(\hbox{\boldmath{\scriptsize$k$}})}{\dot{\bar{A}}_{*}}&(\lambda_{\rm
gw}=+)\\\
\\!\\!\frac{e^{(2)}_{i}\delta\dot{A}_{i*}(\hbox{\boldmath{\scriptsize$k$}})}{\dot{\bar{A}}_{*}}&(\lambda_{\rm
gw}=\times)\end{array}\right.\\!,$

where we have abbreviated the index $\alpha=1$. This indicates that the gauge
field contributes differently to the two linear polarization modes even if
$e_{(1)}^{i}\delta\dot{A}_{i*}$ and $e_{(2)}^{i}\delta\dot{A}_{i*}$ have the
same amplitude [12]. This is because the fluctuation of the amplitude of
$\pi^{i}$, which corresponds to the second term in $\delta\hat{\pi}^{i}$,
includes only the $\lambda=1$ mode under the choice of
$e_{i}^{(2)}\bar{\pi}^{i}_{(1)}=0$ and only the $\lambda=1\,(2)$ mode of the
gauge field contributes to $\lambda_{\rm gw}=+\,(\times)$.

We can similarly compute the sourced gravitational waves for $D^{\prime}\geq
2$ by using Eq. (10). Also for $D^{\prime}\geq 2$ we find that when the
background component of a gauge field is along the direction of $k$
(corresponding to $\Theta_{\alpha}=0$), this gauge field does not contribute
to the gravitational waves, $\gamma^{(\lambda_{\rm gw})}$.

Power spectrum of the gravitational waves.– So far, we have computed the
mapping between the horizon crossing $t_{*}$ and the reheating $t_{f}$ simply
by assuming (A1) and (A2). In what follows, assuming further that the
background anisotropy was still very small at $t=t_{*}$, we compute the power
spectrum of $\delta\gamma_{ij*}$ and $\delta\dot{A}_{i*}$ by adopting the FLRW
background approximation. The amplitudes of the mode functions of the gauge
fields, $\dot{A}^{\lambda}_{(\alpha)}$, are given by

$\displaystyle\frac{|\dot{A}^{\lambda}_{(\alpha)}(t_{*},\,k)|^{2}}{\dot{\bar{A}}_{*(\alpha)}^{2}}=\frac{1}{12}\left(1+\left(\frac{\dot{\bar{f}}_{*}}{H_{*}}\right)^{2}\right)\frac{1}{k^{3}}\frac{\bar{\rho}_{*}}{\bar{\rho}_{A*(\alpha)}}\left(\frac{H_{*}}{M_{\rm
pl}}\right)^{2},$

where $H$ denotes the Hubble parameter. Combining the expressions given above,
we obtain the power spectrum of the primordial gravitational waves for
$D^{\prime}=1$ as

$\displaystyle\langle\gamma^{(\lambda_{\rm
gw})}(t_{f},\,\hbox{\boldmath{$k$}})\gamma^{(\lambda_{\rm
gw})}(t_{f},\,\hbox{\boldmath{$p$}})\rangle$
$\displaystyle=\delta(\hbox{\boldmath{$k$}}+\hbox{\boldmath{$p$}})\frac{2}{k^{3}}\left(\frac{H_{*}}{M_{\rm
pl}}\right)^{2}\\!(1+g_{t}^{\lambda_{\rm gw}}\sin^{2}\Theta)\,.$ (11)

with

$\displaystyle
g_{t}^{\times}=\frac{\bar{\rho}_{*}}{3\bar{\rho}_{A*}}\bar{\hat{{\cal
I}}}_{f}^{2}\left(1+\left(\frac{\dot{\bar{f}}_{*}}{H_{*}}\right)^{2}\right),\quad
g_{t}^{+}=g_{t}^{\times}\cos^{4}\Theta\,.$ (12)

Using $\bar{\hat{{\cal I}}}_{f}\sim 2\Delta\psi\bar{\rho}_{A}^{\rm
max}/\bar{\rho}_{*}$, we find $g_{t}^{\lambda_{\rm
gw}}\sim(\bar{\rho}_{A}^{\rm max}/\bar{\rho}_{*})(\bar{\rho}_{A}^{\rm
max}/\bar{\rho}_{A*})\Delta\psi^{2}$. This indicates that even if the energy
density of the gauge field all the time remains much smaller than the total
energy density, satisfying $(\bar{\rho}_{A}^{\rm
max}/\bar{\rho}_{*})\Delta\psi^{2}\ll 1$, the anisotropic component of the
primordial gravitational waves can be as large as or even larger than the
isotropic one 111In this case, we need to be careful about the terms with
higher powers of $\rho_{A}/\rho$, which are approximately ignored e.g., in
(5).. This is because of the enhancement by $\bar{\rho}_{A}^{\rm
max}/\bar{\rho}_{A*}$, which becomes much larger than 1 when the mode $k$
crosses the horizon before $\rho_{A}$ reaches the maximum value. This is
because the conversion from the fluctuation of the gauge field to the
gravitational waves takes place mainly when $\rho_{A}$ reaches the maximum
value, while the (normalized) power spectrum of the gauge field at the horizon
crossing is inversely proportional to $\bar{\rho}_{A*}$. Because of this
enhancement, even if the gauge field is only sourced by the spectator fields,
which only occupies a small fraction of the total energy density, we can
obtain a large statistical anisotropy [12]. Similarly, we can generalize the
discussion to $D^{\prime}\geq 1$, obtaining

$\displaystyle\langle\gamma^{(\lambda_{\rm
gw})}(t_{f},\,\hbox{\boldmath{$k$}})\gamma^{(\lambda_{\rm
gw})}(t_{f},\,\hbox{\boldmath{$p$}})\rangle$
$\displaystyle=\delta(\hbox{\boldmath{$k$}}+\hbox{\boldmath{$p$}})\frac{2}{k^{3}}\left(\frac{H_{*}}{M_{\rm
pl}}\right)^{2}\\!\left(1+\sum_{\alpha=1}^{D^{\prime}}g_{t(\alpha)}^{\lambda_{\rm
gw}}\sin^{2}\Theta_{\alpha}\right)\,.$ (13)

with $g_{t(\alpha)}^{\lambda_{\rm gw}}\sim(\bar{\rho}_{A(\alpha)}^{\rm
max}/\bar{\rho}_{*})(\bar{\rho}_{A(\alpha)}^{\rm
max}/\bar{\rho}_{A(\alpha)}*)\Delta\psi^{2}$, while the angular dependence of
$g_{t(\alpha)}^{\lambda_{\rm gw}}$ becomes more complicated than the one for
$D^{\prime}=1$.

Adiabatic curvature perturbation.– The (linear) adiabatic perturbation is
given by the summation of the fluctuation of the $e$-folding number and the
longitudinal part of $\delta\gamma_{ij}$ at $t=t_{f}$, where the background
spatial metric is set to $\delta_{ij}$. Under the assumptions (A1) and (A2),
we obtain the dominant contribution of the latter as

$\displaystyle\frac{1}{4}\hat{k}_{i}\hat{k}_{j}\delta\gamma_{ij}(t_{f},\,\hbox{\boldmath{$k$}})=\sum_{\alpha=1}^{D^{\prime}}(\hat{k}_{i}\hat{\bar{\pi}}^{i}_{(\alpha)})^{2}\bar{\hat{{\cal
I}}}_{f(\alpha)}\hat{\bar{\pi}}^{j}_{(\alpha)}\frac{\delta\dot{A}_{j*(\alpha)}(\hbox{\boldmath{$k$}})}{\dot{\bar{A}}_{*}}\,,$

while the former depends significantly on the detail of the models.

For example, when a canonical inflaton field $\phi$ which dominates the energy
density of the Universe does not interact directly with other scalar fields
$\sigma_{I}$ with $I=1,\,\cdots D-1$ nor with the gauge field
$(D^{\prime}=1)$, under the slow-roll approximation, we obtain

$\displaystyle\frac{d\phi}{d\psi}\simeq-M_{\rm
pl}^{2}\frac{V_{\phi}}{V+\rho_{\sigma}+\rho_{A}}\,,$

where $V$ and $V_{\phi}$ denote the scalar potential of $\phi$ and its
derivative and $\rho_{\sigma}$ denotes the energy density of $\sigma_{I}$s.
Assuming that when the total energy density of the gauge field takes the
maximum value $(\psi_{\rm max}\leq\psi\leq\psi_{\rm max}+\Delta\psi)$ both
$\rho_{A}$ and $\rho_{\sigma}$ remain almost constant until $\psi_{f}$, the
above equation can be solved as

$\displaystyle\psi_{f}-\psi_{*}\simeq(\psi_{f}-\psi_{*})_{\phi}+\frac{\rho^{\rm
max}_{\sigma}+\rho^{\rm max}_{A}}{V}\Delta\psi\,,$

where $\rho^{\rm max}_{\sigma}$ denotes $\rho_{\sigma}$ during $\Delta\psi$
but not the maximum value of $\rho_{\sigma}$, and the first term denotes the
$e$-folding number, which is determined only by the dynamics of $\phi$. The
other scalar fields can also provide a angular-independent sub-dominant
contributions to $\psi_{\rm max}-\psi_{*}$, which can be addressed by using
the conventional $\delta N$ formalism. Here and hereafter, we ignore them,
focusing on the leading angular dependent contribution. Perturbing the above
expression, we find that the dominant angular dependent contribution appears
from the fluctuation of $\rho_{\sigma}^{\rm max}$ under the assumptions (A2)
and $\Delta\psi\gg 1$. Expressing the time evolution of $\rho_{A}$ during
$\psi_{*}\leq\psi\leq\psi_{\rm max}$ by using a transfer function ${\cal T}$
as $\rho_{A}^{\rm max}=\cal{T}(\psi_{\rm max}-\psi_{*})\times\rho_{A*}$, we
obtain

$\delta(\psi_{\rm max}-\psi_{*})\simeq-\frac{1}{(\ln{\cal
T})^{\prime}}\frac{\delta\rho_{A*}}{\bar{\rho}_{A*}}\,,$

since we can ignore $\delta\rho_{A}^{\rm max}$ under (A2). Here, the prime
denotes the derivative with respect to the argument of ${\cal T}$, $\psi_{\rm
max}-\psi_{*}$. Using this expression in the fluctuation of
$\rho_{\sigma}^{\rm
max}\simeq\rho_{\sigma*}+(d\rho_{\sigma}/d\psi)_{*}(\psi_{\rm max}-\psi_{*})$,
we obtain

$\displaystyle\zeta(t_{f},\,\hbox{\boldmath{$k$}})$
$\displaystyle\simeq-\frac{V_{*}}{M_{\rm
pl}^{2}V_{\phi*}}\delta\phi_{*}(\hbox{\boldmath{$k$}})-2\frac{\rho_{A}^{\rm
max}}{\rho_{*}}\Delta\psi\frac{\hat{\bar{\pi}}^{i}\delta\dot{A}_{i*}(\hbox{\boldmath{$k$}})}{\dot{\bar{A}}_{*}}$
$\displaystyle\qquad\qquad\qquad\quad\times\left(\frac{d\rho_{\sigma}}{d\psi}\\!\\!\Big{|}_{*}\frac{1}{(\ln{\cal
T})^{\prime}\rho_{A}^{\rm max}}+\cos^{2}\Theta\right)\\!.$

The second term in the parenthesis is the model independent contribution which
comes from the longitudinal mode of $\delta\gamma_{ij}$ and the first one
comes from the fluctuation of the $e$-folding number. The different
$\rho_{A*}$ results in the difference in $\rho_{\sigma}^{\rm max}$, since the
$e$-folding number spent until $\psi_{\rm max}$ differs, generating the
fluctuation of the $e$-folding number, $\delta(\psi_{f}-\psi_{*})$. For the
specific choice of $f$ and the potential of $\sigma_{I}$ with $D=2$, the above
formula of $\zeta$ reproduces the result obtained in Ref. [12]. Since the
statistical anisotropy in the power spectrum of $\zeta$ is suppressed by the
slow-roll parameter $\varepsilon=(M_{\rm pl}V_{\phi}/V)^{2}/2$ compared to
that of the gravitational waves $g_{t}^{\lambda_{\rm gw}}$, the statistical
anisotropy in the power spectrum of the gravitational waves can be as large as
${\cal O}(1)$, keeping the scalar spectrum consistent with the present
observations [12]. This remains the same as long as the amplitude of the first
term in the parenthesis of $\zeta$ does not exceed ${\cal O}(1)$. Using these
formulae, we can also compute the cross-correlation of $\zeta$ and
$\gamma^{+}$, which takes a non-vanishing value.

Primordial non-Gaussianity.– Using the g$\delta N$ formalism, the local type
non-Gaussianity can be computed easily. Here, let us only provide the order
estimation, leaving a detailed computation for elsewhere. We obtain the local
type non-Gaussianity of $\gamma^{(\lambda_{\rm gw})}$ sourced by the gauge
fields as

$\displaystyle f_{\rm
NL}^{\gamma,A}\sim\sum_{\alpha=1}^{D^{\prime}}\Bigl{(}\frac{\bar{\rho}_{*}}{\bar{\rho}_{A(\alpha)*}}\Bigr{)}^{\\!\\!2}\,\bar{\hat{{\cal
I}}}_{f(\alpha)}^{3}\sim\sum_{\alpha=1}^{D^{\prime}}\Bigl{(}\frac{\bar{\rho}_{A(\alpha)}^{\rm
max}}{\bar{\rho}_{A(\alpha)*}}\Bigr{)}^{\\!\\!2}\frac{\bar{\rho}_{A(\alpha)}^{\rm
max}}{\bar{\rho}_{*}}\Delta\psi_{(\alpha)}^{3}\,,$

which can be also enhanced by the square of $\bar{\rho}_{A(\alpha)}^{\rm
max}/\bar{\rho}_{A(\alpha)*}$, which can be very large when the corresponding
gauge field reaches the maximum value after the horizon crossing. Here, we
have ignored the angular dependence which requires a more detailed
computation. The local type non-Gaussianity of $\zeta$ which stems from the
longitudinal part of $\delta\gamma_{ij}$ is suppressed by $\varepsilon^{2}$
compared to $f_{\rm NL}^{\gamma,A}$, while there is also model dependent
contribution in the fluctuation of the $e$-folding number.

Summary.– In this letter, we showed the g$\delta N$ formalism can largely
facilitate the computation of $\zeta$ and the gravitational waves by
considering an inflation model with U(1) gauge fields. Since the g$\delta N$
formalism generically applies to a model with the locality and the spatial
diffeomorphism invariance, various applications will be possible.

###### Acknowledgements.

T. T. is supported by Grant-in-Aid for Scientific Research (A) and (C) under
Contract Nos. JP23H00110 and JP20K03928, respectively. Y. U. is supported by
Grant-in-Aid for Scientific Research (B) under Contract Nos. JP19H01894 and
JP23H01177 and Fostering Joint International Research (B) under Contract No.
21KK0050.

## References

  * Starobinsky [1982] A. A. Starobinsky, Dynamics of Phase Transition in the New Inflationary Universe Scenario and Generation of Perturbations, Phys. Lett. 117B, 175 (1982).
  * Sasaki and Stewart [1996] M. Sasaki and E. D. Stewart, A General analytic formula for the spectral index of the density perturbations produced during inflation, Prog. Theor. Phys. 95, 71 (1996), arXiv:astro-ph/9507001 [astro-ph] .
  * Sasaki and Tanaka [1998] M. Sasaki and T. Tanaka, Superhorizon scale dynamics of multiscalar inflation, Prog. Theor. Phys. 99, 763 (1998), arXiv:gr-qc/9801017 [gr-qc] .
  * Salopek and Bond [1990] D. S. Salopek and J. R. Bond, Nonlinear evolution of long wavelength metric fluctuations in inflationary models, Phys. Rev. D42, 3936 (1990).
  * Wands  _et al._ [2000] D. Wands, K. A. Malik, D. H. Lyth, and A. R. Liddle, A New approach to the evolution of cosmological perturbations on large scales, Phys. Rev. D62, 043527 (2000), arXiv:astro-ph/0003278 [astro-ph] .
  * Tanaka and Urakawa [2021] T. Tanaka and Y. Urakawa, Anisotropic separate universe and Weinberg’s adiabatic mode, JCAP 07, 051, arXiv:2101.05707 [astro-ph.CO] .
  * Shibata and Sasaki [1999] M. Shibata and M. Sasaki, Black hole formation in the Friedmann universe: Formulation and computation in numerical relativity, Phys. Rev. D60, 084002 (1999), arXiv:gr-qc/9905064 [gr-qc] .
  * [8] T. Tanaka and Y. Urakawa, in preparation, .
  * Watanabe  _et al._ [2009] M.-a. Watanabe, S. Kanno, and J. Soda, Inflationary Universe with Anisotropic Hair, Phys. Rev. Lett. 102, 191302 (2009), arXiv:0902.2833 [hep-th] .
  * Watanabe  _et al._ [2010] M.-a. Watanabe, S. Kanno, and J. Soda, The Nature of Primordial Fluctuations from Anisotropic Inflation, Prog. Theor. Phys. 123, 1041 (2010), arXiv:1003.0056 [astro-ph.CO] .
  * Fujita  _et al._ [2019] T. Fujita, R. Namba, and I. Obata, Mixed Non-Gaussianity from Axion-Gauge Field Dynamics, JCAP 04, 044, arXiv:1811.12371 [astro-ph.CO] .
  * Fujita  _et al._ [2018] T. Fujita, I. Obata, T. Tanaka, and S. Yokoyama, Statistically Anisotropic Tensor Modes from Inflation, JCAP 07, 023, arXiv:1801.02778 [astro-ph.CO] .
  * Note [1] In this case, we need to be careful about the terms with higher powers of $\rho_{A}/\rho$, which are approximately ignored e.g., in (5).