Copyright for this paper by its authors. Use permitted under Creative Commons
License Attribution 4.0 International (CC BY 4.0).

ISWC 2023: Deep Learning for Knowledge Graph(DL4KG), November 06-10, 2023,
Athens, Greece

[orcid=0000-0002-1851-5390<EMAIL_ADDRESS>] [1]

[orcid=0009-0004-6298-5711<EMAIL_ADDRESS>]

[1]Corresponding author.

# NNKGC: Improving Knowledge Graph Completion with Node Neighborhoods

Irene Li  The University of Tokyo, Tokyo, Japan  Boming Yang

(2023)

###### Abstract

Knowledge graph completion (KGC) aims to discover missing relations of query
entities. Current text-based models utilize the entity name and description to
infer the tail entity given the head entity and a particular relation.
Existing approaches also consider the neighborhood of the head entity.
However, these methods tend to model the neighborhood using a flat structure
and are only restricted to 1-hop neighbors. In this work, we propose a node
neighborhood-enhanced framework for knowledge graph completion. It models the
head entity neighborhood from multiple hops using graph neural networks to
enrich the head node information. Moreover, we introduce an additional edge
link prediction task to improve KGC. Evaluation on two public datasets shows
that this framework is simple yet effective. The case study also shows that
the model is able to predict explainable predictions.

###### keywords:

Knowledge Graph BERT Link Prediction Graph Convolutional Networks

##  1 Introduction

A knowledge graph (KG) is a structured representation of entities, their
properties, and their relations [1]. It has become increasingly popular in
recent years due to its ability to support various applications, including
question answering, recommendation systems, and semantic search [2, 3, 4].
However, a KG is usually defined by domain experts, so a few relations might
be missing. Knowledge Graph Completion (KGC) aims to discover these missing
relations in a given KG [5].

The methods for KGC mainly fall into two groups: embedding-based and text-
based. Embedding-based methods consider the knowledge graph relations as graph
structures, then learn low-rank embeddings to represent entities and
relations, such as TransE [6], TuckER [7] and so on. However, such methods
fail in any inductive scenarios. Text-based methods infer relations based on
entities and the corresponding descriptions through representation learning to
solve this problem [8, 9, 10]. KEPLER [11] utilizes pretrained language
modeling representation with knowledge embedding to integrate factual
knowledge for KGC. The recent SimKGC [12] model applies a contrastive learning
approach with new methods for negative sampling, outperforming several
embedding-based methods.

Most existing work typically models the head entity $h$ and the query relation
$r$ together, then focuses on the modeling for predicting the relation for the
tuple $(h,r,?)$ to find the correct tail entity $t$. The modeling for $(h,r)$
tends only to consider the head entity name or descriptions and the query
relation name. Besides, existing approaches also consider the neighborhood of
the head entity [13, 14]. However, these methods tend to model the
neighborhood using a flat structure so that a simple attention-based method
can be applied; or they are only restricted to 1-hop neighbors for
computational consideration. In this work, we propose considering the
neighborhood from multiple hops to improve the richness of the head entity
information for further KGC predictions. By doing so, the head entity contains
much explicit information that may help the relation prediction, especially
when the head entity name/description is short and not informative. We model
the neighborhood using a graph neural network, making it easy to propagate
information within multiple hops away. Besides, to further enhance the
correlation with these neighbors and the tail entity, we add an extra link
prediction loss to the KGC task, and we show this could improve the KGC
performance.

Our main contribution of this work is to propose a node neighborhood-enhanced
model for knowledge graph completion by modeling the neighbors with graph
neural networks [15]. Moreover, the experiment shows that this framework is
effective in two public datasets. Additional case studies show that this model
could provide explainable predictions. We release our code at
https://github.com/IreneZihuiLi/NNKGC.

##  2 Method

We propose the Node Neighborhood-enhanced model for knowledge Graph Completion
(NNKGC), which consists of the node encoder and the relation prediction
module, as shown in Fig. 1. Following the setting of the previous work [12,
8], we define the task to be the following: given an incomplete knowledge
graph $\mathcal{G}$, we need to find the tail entity $t$ by providing the head
entity $h$ and the query relation $r$: $(h,r,?)$.

Node Encoder Typically, we model the head and tail entity nodes separately. In
our case, each entity node contains a node name in a short phrase or a keyword
(name), and a node description in a sentence (desc). The relation $r$ is
usually represented by a phrase or a keyword (rel). To incorporate the text
features into the knowledge graph, we utilize the pre-trained language models
(PLMs), BERT [16], to encode the head entity and relations. Each entity node
is represented by:

$e_{node}=\mathbf{BERT}[name,desc,\texttt{[SEP]},rel]$ (1)

To further consider the head node neighborhood to enrich the encoding, we
apply a graph neural network [17] to model the neighbors. We first collect all
the neighbors within $k$-hop as a list of nodes: $nb1,nb2,...nbl$. Then we
encode the head entity and relation as follows:

$\displaystyle X$ $\displaystyle=[e_{head},e_{nb1},e_{nb2},...e_{nbl}],$ (2)
$\displaystyle e_{hr}$ $\displaystyle=\mathbf{GNN}(X,A),$ (3)

In which, $X$ denotes the node embeddings of the entities within the
neighborhood, including the head node, as well as the neighbors, and $A$ is
the adjacency matrix indicating the directed connections among the entities.
We define the neighborhood of a head node to be the connected entities,
including incoming and outgoing connections. In Fig 1, we illustrate 1-hop
(blue) and 2-hop (grey) neighbors and ignore outgoing neighbors for
simplicity. There is a number of choices for the graph encoder, we investigate
three types: graph convolutional network (GCN) [18], graph attention network
(GAT) [19], and GraphSAGE [20].

Figure 1: NNKGC Model Illustration: we initially identify the neighborhood of
the head node. Following this, we employ a graph encoder to encode both the
head node and its corresponding relation. Simultaneously, the tail node is
encoded. For the final link prediction, we utilize cosine similarity. A novel
aspect of our methodology is the introduction of the VGAE loss. This is
optimized in conjunction with the knowledge graph task loss during the model’s
training phase.

Relation Prediction To make the prediction on the tail entities, we followed
the recent SimKGC [12] model. SimKGC introduces a contrastive learning
approach to improve the negative sampling and applies the InfoNCE loss as the
loss function. The prediction is simply a score based on cosine similarity
between $e_{hr}$ and $e_{t}$. Then the tail entities are ranked by this score
and the largest one is chosen as the prediction. The loss for the tail entity
prediction is denoted as ${L}_{kg}$.

Neighborhood Edge Prediction To further strengthen the effect of neighborhood
entities for the KGC task, we ask the model to predict neighborhood edges
(i.e., the highlighted edge in the figure). In other words, we randomly mask
some edges in matrix $A$, and let the model reconstruct these masked
relations. There are multiple ways to achieve this, and we apply a variational
graph autoencoder (VGAE) [21] to conduct the missing link prediction. It is
shown to be effective for several NLP tasks [21, 22]. Specifically, we apply
vanilla GCNs to model the mean and standard deviation of the node features,
then conduct dot-product to predict the existence of a given node pair. Note
that this is only done for the neighborhood. It will generate a loss on the
reconstructed edges $\mathcal{L}_{edge}$. So the final loss becomes:

$\mathcal{L}=\lambda\mathcal{L}_{kg}+(1-\lambda)\mathcal{L}_{edge}$ (4)

Here, $\lambda$ is an empirical constant to control the effect of the edge
prediction, and it ranges between 0 and 1.

Table 1: Dataset statistics: vocab indicates the vocabulary size; node mean
indicates the average number of tokens for the node descriptions, and the node
median indicates the median number of tokens.

dataset entity relation mean median train/valid/test vocab node mean node
median WN18RR 40,943 11 3.7 2 86, 835/3,034/3,134 18,741 17.1 15 FB15k-237
14,541 237 2.8 2 272,115/17,535/20,466 28,538 114.6 76

##  3 Experiment

Shown in Tab. 1, we conduct knowledge graph completion on two public datasets:
FB15k-237 [23], and WN18RR [24]. We follow the previous work [11, 12] to
report the following evaluation metrics: mean reciprocal rank (MRR),
Hits@1,3,10.

Hyperparameters The BERT encoder is the pretrained bert-base-uncased model
111https://huggingface.co/bert-base-uncased. AdamW optimizer [25] was applied
for training. For modeling the neighborhood, we typically apply two layers of
neural networks, and the dimension for the graph node representation is 768.
We set the number of attention heads for the GAT encoder to 3. We set the
$\lambda$ to 0.2. We conduct experiments on 4 A100 GPUs (with 40GB memory).
Training on WN18RR with 30 epochs took about 3 hours; training on FB15k237
with 5 epochs took less than one hour. More experimental hyperparameters can
be found in the code URL.

Table 2: Main results on selected datasets. We report MRR, H@1, H@3 and H@10
scores.

FB15k-237 | MRR | H@1 | H@3 | H@10
---|---|---|---|---
TransE [6] | 27.9 | 19.8 | 37.6 | 44.1
TuckER [7] | 35.8 | 26.6 | 39.4 | 54.4
MTL-KGC [26] | 26.7 | 17.2 | 29.8 | 45.8
SimKGC [12] | 33.6 | 24.9 | 36.2 | 51.1
1-5 NNKGC | 33.2 | 24.6 | 35.9 | 50.4
NNKGCe | 33.8 | 25.2 | 36.5 | 51.5
WN18RR | MRR | H@1 | H@3 | H@10
TransE [6] | 24.3 | 4.3 | 44.1 | 53.2
TuckER [7] | 47.0 | 44.3 | 48.2 | 52.6
MTL-KGC [26] | 33.1 | 20.3 | 38.3 | 59.7
SimKGC [12] | 66.6 | 58.7 | 71.7 | 80.0
1-5 NNKGC | 67.3 | 59.8 | 71.7 | 80.5
NNKGCe | 67.4 | 59.6 | 72.2 | 81.2

Main Results We compare with two embedding-based methods: TransE [6], and
TuckER [7]. TransE learns low-dimensional embeddings of the entities, and
TuckER applies Tucker decomposition of the binary tensor representation for
knowledge graph triples. We also include text-based methods for baselines.
MTL-KGC [26] is a multi-task learning method for KGC, and SimKGC [12]. We
compare two settings of our NNKGC model: without the edge loss
$\mathcal{L}_{edge}$ (NNKGC) and with the edge loss (NNKGCe). Text-based
methods performs better than embedding-based methods. Without the VGAE loss,
our NNKGC is competitive with the best baseline, SimKGC. In general, our best
model NNKGCe outperforms the selected baselines in both datasets.

##  4 Ablation Study

Neighborhood Encoder We then study different neighborhood graph encoders in
Eq. 3. We compare three widely-used graph models: GCN, GAT, and GraphSAGE. In
Fig. 2, we report the average performance on the two datasets. We could
observe that GraphSAGE and GAT are slightly better than the vanilla GCN
encoder.

WNFB$0.2$$0.4$$0.6$$0.8$Avg.GCNGraphSAGEGAT Figure 2: Comparison of multiple
neighborhood graph modeling: GCN, GraphSAGE and GAT.

We provide a detailed evaluation of various graph encoders for modeling the
head node neighborhood, shown in Tab. 4 and 3. In general, the GCN encoder is
worse than the other two. In general, GraphSAGE shows better stability in both
datasets.

Table 3: Comparison of multiple neighborhood graph modeling: GCN, GraphSAGE, and GAT on the FB15k-237 dataset. Encoder | MRR | H@1 | H@3 | H@10
---|---|---|---|---
GCN | 32.7 | 23.8 | 35.4 | 50.6
GraphSAGE | 33.8 | 25.2 | 36.5 | 51.5
GAT | 33.7 | 25.1 | 36.2 | 51.2
Table 4: Comparison of multiple neighborhood graph modeling: GCN, GraphSAGE, and GAT on the WN18RR dataset. Encoder | MRR | H@1 | H@3 | H@10
---|---|---|---|---
GCN | 65.4 | 57.0 | 70.7 | 80.2
GraphSAGE | 65.8 | 56.4 | 72.4 | 81.7
GAT | 67.4 | 59.6 | 72.2 | 81.2
Table 5: An example from FB15k-237: we show the ground truth tail and
predictions, as well as a selection of neighbor entities. We highlight direct
triggering tokens and other relevant tokens.

Head: Star Wars Episode IV: A New Hope
---
Head node description: Star Wars, later retitled Star Wars Episode IV: A New
Hope, is a 1977 American epic space opera film written and directed by George
Lucas…
Relation: nominated for (an award)
Ground Truth Tail: Academy Award for Best Sound Mixing
Predicted Tails: Academy Award for Best Sound Mixing (0.71), BAFTA Award for
Best Special Visual Effects (0.656), Golden Globe Award for Best Original
Score (0.65)
Neighbor 1, Name: Academy Award for Best Production Design
Description: The Academy Awards are the oldest awards ceremony for
achievements in motion pictures. The Academy Award for Best Production Design
recognizes achievement in art direction on a film. The category’s original
name was Best Art Direction, but was changed to its current name in 2012 for
the 85th Academy Awards.
Node relation: nominated for (an award)
1-1 Neighbor 2, Name: Ben Burtt
Description: Benjamin Ben Burtt, Jr. is an American sound designer, film
editor, director, screenwriter, and voice actor. He has worked as sound
designer on various films including: the Star Wars and Indiana Jones film
series…
Node relation: winner, won (an award)
1-1 Neighbor 3, Name: John Williams
Description: John Towner Williams is an American composer, conductor and
pianist…
Node relation: music film

Neighborhood Hops As we focus on the neighborhood, we now study the effect of
the neighborhood scale. We conduct experiments on 1-hop, 2-hop and 3-top
neighbors to the head entity. Shown in Fig. 3, we compare the evaluation
metrics as well as the average performance (Avg. on the WN18RR dataset.

Figure 3: The effect of hop numbers, results on WN18RR. $1$$2$$3$$0.6$$0.8$Hop
NumberMRRH@1H@3H@10Avg.

We also show the bar chart of how the number of hops affects the performance
for FB15k-237 in Fig. 4. One may observe a similar trend to that of WN18RR. As
more hops will bring more neighboring entities, some noises may be added,
which may make the performance worse.

Figure 4: Effect of the Number of Hops, results on
FB15k-237.$1$$2$$3$$0.2$$0.3$$0.4$Number of HopMRRH@1H@3H@10AVG

In general, integrating more neighbors is not contributing to the performance.
Especially, H@1 drops significantly when the hop number increases. A possible
reason might be that more neighbors may introduce more noises which confuses
the knowledge graph completion task. Moreover, doing so brings more
computational burden to the model as more nodes participate in the graph
modeling part. We did not conduct experiments on 4 and more hops, as the total
training time is incredibly long.

Table 6: An example from WN18RR: we show the ground truth tail and
predictions, as well as a selection of neighbor entities. We highlight direct
triggering tokens and other relevant tokens.

Head: grant (NN)
---
Head node description: any monetary aid
Relation: hypernym
Ground Truth Tail: financial aid (NN1)
Predicted Tails: financial aid (NN1) (0.81), grant (NN1) (0.68), foreign aid
(NN1)(0.572),
Neighbor 1, Name: grant in aid (NN2)
Description: a grant to a person or school for some educational project
Node relation: hypernym
1-1 Neighbor 2, Name: grant (VB2)
Description: give as judged due or on the basis of merit; the referee awarded
a free kick to the team; the jury awarded a million dollars to the plaintiff;
funds are granted to qualified researchers.
Node relation: derivationally related form

Table 7: An example from FB15k-237: we show the ground truth tail and
predictions, as well as a selection of neighbor entities. We highlight direct
triggering tokens and other relevant tokens.

Head: Jack London
---
Head node description: John Griffith Jack London was an American author,
journalist, and social activist. He was a pioneer in the then-burgeoning world
of commercial magazine fiction and was one of the first fiction writers to
obtain worldwide celebrity and a large fortune from his fiction alone…
Relation: nationality person people
Ground Truth Tail: United States of America
Predicted Tails: United States of America (0.91), Confederate States of
America (0.40), Union (0.34),
Neighbor 1, Name: San Francisco
Description: San Francisco, officially the City and County of San Francisco,
is the leading financial and cultural center of Northern California and the
San Francisco Bay Area…
Node relation: place of birth person people
1-1 Neighbor 2, Name: Oakland
Description: Oakland, located in the U.S. state of California, is a major West
Coast port city and the busiest port for San Francisco Bay and all of Northern
California. It is the third largest city in the San Francisco Bay Area, the
eighth-largest city in the state, and the 47th-largest city
Node relation: location place, lived places
1-1 Neighbor 3, Name: University of California, Berkeley
Description: The University of California, Berkeley, is a public research
university located in Berkeley, California, United States. The university
occupies 1,232 acres on the eastern side of the San Francisco Bay with the
central campus resting on 178 acres. Berkeley is the flagship institution of
the 10 campus University of California
Node relation:students, graduates educational institution

Case Study We select an example from FB15k-237 in Tab. 4, and study how the
neighbor entities and corresponding relations help with the KGC task. Given
the head entity to be the movie Star Wars IV, and the relation to be nominated
for (an award), our model predicts a ranked possible tail entity list. The
entity Academy Award for Best Sound Mixing achieves the highest score, which
is precisely the ground truth. We also randomly select three 1-hop neighbors,
and their relations with the head entity in the training set. As we can see,
the first neighbor is one of the Academy Award categories, and the relation is
the same as the query relation. The second neighbor is a sound designer. The
description shows that this person worked on the Star Wars series. The third
neighbor is a composer who participates in the music production for the movie.
We highlight some potential triggering tokens in the table. As we can see,
some tokens are highly related to the ground truth tail entity, i.e., neighbor
1, the name  Academy Award. Besides, the token composer and music are also
semantically related to the ground truth tail entity name Best Sound Mixing.
From this case study, one can notice that the neighborhood has some positive
effects on the KGC task, which is consistent with our motivation. Moreover, we
show that our framework has the potential to explain the predictions, which is
an essential step for further fact verification.

We present two more case studies. Tab. 4 shows a random example from WN18RR.
As we can see, the entity description tends to be shorter, but we can still
find text chunks that are highly related to the ground truth tail entity
(financial aid and funds). Similarly, Tab. 4 gives an example of finding a
person’s nationality. We can observe that the neighbors are some relevant
locations about this head entity person (i.e., lives in Oakland, U.S. state of
California), which contributes to the final prediction on the tail entity
(United States of America).

##  5 Conclusion

In this work, we proposed a node neighborhood-enhanced model for knowledge
graph completion, by modeling the neighbors with graph neural networks. We
showed that the framework is simple but effective. Case studies also show that
it is possible to provide explainable results.

## References

  * Ehrlinger and Wöß [2016] L. Ehrlinger, W. Wöß,  Towards a definition of knowledge graphs.,  SEMANTiCS (Posters, Demos, SuCCESS) 48 (2016) 2.
  * Steiner et al. [2012] T. Steiner, R. Verborgh, R. Troncy, J. Gabarro, R. Van de Walle,  Adding realtime coverage to the google knowledge graph,  in: 11th International Semantic Web Conference (ISWC 2012), volume 914, Citeseer, 2012, pp. 65–68.
  * Wang et al. [2021] X. Wang, T. Huang, D. Wang, Y. Yuan, Z. Liu, X. He, T.-S. Chua,  Learning intents behind interactions with knowledge graph for recommendation,  Proceedings of the Web Conference 2021 (2021).
  * Liang et al. [2022] L. Liang, Y. Li, M. Wen, Y. Liu,  Kg4py: A toolkit for generating python knowledge graph and code semantic search,  Connection Science 34 (2022) 1384 – 1400.
  * Lin et al. [2015] Y. Lin, Z. Liu, M. Sun, Y. Liu, X. Zhu,  Learning entity and relation embeddings for knowledge graph completion,  in: AAAI Conference on Artificial Intelligence, 2015\.
  * Bordes et al. [2013] A. Bordes, N. Usunier, A. García-Durán, J. Weston, O. Yakhnenko,  Translating embeddings for modeling multi-relational data,  in: NIPS, 2013.
  * Balazevic et al. [2019] I. Balazevic, C. Allen, T. Hospedales,  TuckER: Tensor factorization for knowledge graph completion,  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), Association for Computational Linguistics, Hong Kong, China, 2019, pp. 5185–5194. URL: https://aclanthology.org/D19-1522. doi:10.18653/v1/D19-1522.
  * Yao et al. [2019] L. Yao, C. Mao, Y. Luo,  Kg-bert: Bert for knowledge graph completion,  ArXiv abs/1909.03193 (2019).
  * Kim et al. [2020] B. Kim, T. Hong, Y. Ko, J. Seo,  Multi-task learning for knowledge graph completion with pre-trained language models,  in: Proceedings of the 28th International Conference on Computational Linguistics, International Committee on Computational Linguistics, Barcelona, Spain (Online), 2020, pp. 1737–1743. URL: https://aclanthology.org/2020.coling-main.153. doi:10.18653/v1/2020.coling-main.153.
  * Safavi and Koutra [2021] T. Safavi, D. Koutra,  Relational world knowledge representation in contextual language models: A review,  ArXiv abs/2104.05837 (2021).
  * Wang et al. [2021] X. Wang, T. Gao, Z. Zhu, Z. Zhang, Z. Liu, J. Li, J. Tang,  KEPLER: A unified model for knowledge embedding and pre-trained language representation,  Transactions of the Association for Computational Linguistics 9 (2021) 176–194. URL: https://aclanthology.org/2021.tacl-1.11. doi:10.1162/tacl_a_00360.
  * Wang et al. [2022] L. Wang, W. Zhao, Z. Wei, J. Liu,  Simkgc: Simple contrastive knowledge graph completion with pre-trained language models,  ArXiv (2022).
  * Zhuo et al. [2022] J. Zhuo, Q. Zhu, Y. Yue, Y. Zhao, W. Han,  A neighborhood-attention fine-grained entity typing for knowledge graph completion,  Proceedings of the Fifteenth ACM International Conference on Web Search and Data Mining (2022).
  * Borrego et al. [2021] A. Borrego, D. Ayala, I. Hernández, C. R. Rivero, D. Ruiz,  Cafe: Knowledge graph completion using neighborhood-aware features,  Eng. Appl. Artif. Intell. 103 (2021) 104302.
  * Wu et al. [2019] Z. Wu, S. Pan, F. Chen, G. Long, C. Zhang, P. S. Yu,  A comprehensive survey on graph neural networks,  IEEE Transactions on Neural Networks and Learning Systems 32 (2019) 4–24.
  * Devlin et al. [2019] J. Devlin, M.-W. Chang, K. Lee, K. Toutanova,  BERT: Pre-training of deep bidirectional transformers for language understanding,  in: Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers), Association for Computational Linguistics, Minneapolis, Minnesota, 2019, pp. 4171–4186. URL: https://aclanthology.org/N19-1423. doi:10.18653/v1/N19-1423.
  * Scarselli et al. [2009] F. Scarselli, M. Gori, A. C. Tsoi, M. Hagenbuchner, G. Monfardini,  The graph neural network model,  in: IEEE International Joint Conference on Neural Networks, IEEE, 2009, pp. 2187–2192.
  * Kipf and Welling [2016] T. Kipf, M. Welling,  Semi-supervised classification with graph convolutional networks,  ArXiv abs/1609.02907 (2016).
  * Velickovic et al. [2017] P. Velickovic, G. Cucurull, A. Casanova, A. Romero, P. Lio’, Y. Bengio,  Graph attention networks,  ArXiv abs/1710.10903 (2017).
  * Hamilton et al. [2017] W. L. Hamilton, Z. Ying, J. Leskovec,  Inductive representation learning on large graphs,  in: NIPS, 2017.
  * Li et al. [2021] I. Li, V. Yan, T. Li, R. Qu, D. Radev,  Unsupervised cross-domain prerequisite chain learning using variational graph autoencoders,  in: Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing (Volume 2: Short Papers), 2021, pp. 1005–1011.
  * Xie et al. [2021] Q. Xie, J. Huang, P. Du, M. Peng, J. Nie,  Inductive topic variational graph auto-encoder for text classification,  in: North American Chapter of the Association for Computational Linguistics, 2021.
  * Toutanova et al. [2015] K. Toutanova, D. Chen, P. Pantel, H. Poon, P. Choudhury, M. Gamon,  Representing text for joint embedding of text and knowledge bases,  in: Conference on Empirical Methods in Natural Language Processing, 2015.
  * Dettmers et al. [2017] T. Dettmers, P. Minervini, P. Stenetorp, S. Riedel,  Convolutional 2d knowledge graph embeddings,  in: AAAI Conference on Artificial Intelligence, 2017\.
  * Loshchilov and Hutter [2017] I. Loshchilov, F. Hutter,  Decoupled weight decay regularization,  in: International Conference on Learning Representations, 2017.
  * Kim et al. [2020] B. Kim, T. Hong, Y. Ko, J. Seo,  Multi-task learning for knowledge graph completion with pre-trained language models,  in: International Conference on Computational Linguistics, 2020.


# Learning deep autoregressive models for hierarchical data

Carl R. Andersson  Niklas Wahlström  Thomas B. Schön  Department of
information technology, Uppsala University, Sweden, (email:
<EMAIL_ADDRESS> Carl R. Andersson
Niklas Wahlström  Thomas B. Schön  Department of information technology,
Uppsala University, Sweden, (email:
<EMAIL_ADDRESS>
# Learning deep autoregressive models for hierarchical data

Carl R. Andersson  Niklas Wahlström  Thomas B. Schön  Department of
information technology, Uppsala University, Sweden, (email:
<EMAIL_ADDRESS> Carl R. Andersson
Niklas Wahlström  Thomas B. Schön  Department of information technology,
Uppsala University, Sweden, (email:
<EMAIL_ADDRESS>
###### Abstract

We propose a model for hierarchical structured data as an extension to the
stochastic temporal convolutional network. The proposed model combines an
autoregressive model with a hierarchical variational autoencoder and
downsampling to achieve superior computational complexity. We evaluate the
proposed model on two different types of sequential data: speech and
handwritten text. The results are promising with the proposed model achieving
state-of-the-art performance.

###### Abstract

We propose a model for hierarchical structured data as an extension to the
stochastic temporal convolutional network. The proposed model combines an
autoregressive model with a hierarchical variational autoencoder and
downsampling to achieve superior computational complexity. We evaluate the
proposed model on two different types of sequential data: speech and
handwritten text. The results are promising with the proposed model achieving
state-of-the-art performance.

###### keywords:

Deep learning, variational autoencoders, nonlinear systems

Please cite this version:

Carl R. Andersson, Niklas Wahlström, Thomas B. Schön. Learning deep
autoregressive models for hierarchical data. In Proceedings of 19th IFAC
Symposium on System Identification (SYSID), Padova, Italy (online), 2021.

⬇ @inproceedings{CarlAndersson2021, author = {Carl R. Andersson and Niklas
Wahlstr\"{o}m and Thomas Bo Sch\"{o}n}, title = {Learning deep autoregressive
models for hierarchical data}, booktitle = {Proceedings of 19th IFAC Symposium
on System Identification (SYSID)}, address = {Padova, Italy (online)}, year =
{2021} }

††thanks:  This research was financially supported by the projects _Learning
flexible models for nonlinear dynamics_ (contract number: 2017-03807) by the
Swedish Research Council, by _AI4Research_ at Uppsala University and by _Kjell
och Märta Beijer Foundation_.

##  1 Introduction

System identification and sequence modeling with deep learning are two
different research areas that essentially solve the same problem, to model a
sequence, $y_{1:T}$, for predictive purposes. Even though the methods are very
similar there are differences when it comes to what kind of data the method
typically is applied to. Whereas system identification traditionally has
focused on systems with relatively short memory and small datasets (e.g.
systems on the nonlinear benchmark website (Schoukens, 2020)), deep learning
has focused on the opposite, i.e. systems with long memory and larger datasets
(e.g. text modeling). This dichotomy is not a product of any fundamental
difference in methodology, but rather an effect of the rapid advancement in
computational power that has accompanied the advancements in deep learning.

Neural networks are by no means new to the system identification community,
they have in fact a rather rich history (Sjöberg et al., 1995). What deep
learning brought with it was instead a new paradigm of model design and
implicit regularization. Whereas the system identification community typically
view neural networks as black box function approximators, the deep learning
community instead augment how these black box models are built up and adapt
them to the data, e.g. long-short term memory (LSTM) models (Hochreiter and
Urgen Schmidhuber, 1997), convolutional neural networks (CNN) (Krizhevsky et
al., 2012) and temporal convolutional networks (TCN). A very prominent example
of the connection between model and data structure is indeed deep CNNs applied
to images, where the model mirrors the hierarchical composition of natural
images and locality of low-level features (LeCun et al., 2015) while high-
level features are represented in a downsampled image. Such hierarchical
composition can also be seen in many sequential datasets, e.g. language
modeling, handwritten text modeling (see Figure 1(a)), human motion tracking
or other systems of (hierarchical) switching nature where the high-level
features evolve at a slower pace than low-level features (Chung et al., 2016;
Koutník et al., 2014).

(a) Dataset samples

(b) Generated samples

Figure 1: An illustration of samples from the dataset and generated samples
from our model for the IAM-OnDB dataset (Liwicki and Bunke, 2005).

In the light of this we propose a sequential hierarchical latent variable
model, represented in Figure 2(a), that focus on modeling multiscale
hierarchically structured sequential data. This can seen as an extension the
Stochastic TCN (Aksan and Hilliges, 2019) which also uses a sequential
hierarchical latent variable model although without the multiscale property.
The latent variables $z=\\{z^{(l)}_{1:T}\\}$ are structured sequentially with
hierarchical layers, where layer $l$ is downsampled with a factor $S^{(l)}$.
To facilitate training we propose to relax the state space formulation
similarily to Aksan and Hilliges (2019), see Figure 2(b). With this we get a
deep sequential model that exploits the hierarchical nature of the data
similarly to how a deep convolutional network exploits the nature of images.
We realize the model using a combination of a Wavenet model (Van Den Oord et
al., 2016) for the sequential component with a hierarchical variational
autoencoder (Rezende et al., 2014; Kingma and Welling, 2014) for the latent
variables. We employ the proposed model on two different datasets and show
that by including the prior information of the hierarchical nature of the data
we can improve the model in terms of both parameter effectiveness and
computational complexity.

pt$z^{(3)}_{t}$$z^{(2)}_{t}$$z^{(2)}_{t+2}$$z^{(1)}_{t}$$z^{(1)}_{t+1}$$z^{(1)}_{t+2}$$z^{(1)}_{t+3}$$y_{t}$$y_{t+1}$$y_{t+2}$$y_{t+3}$$\dots$$\dots$pt$\dots$pt$\dots$pt$\dots$pt$\dots$

(a) State space form

pt$z^{(3)}_{t}$$z^{(2)}_{t}$$z^{(2)}_{t+2}$$z^{(1)}_{t}$$z^{(1)}_{t+1}$$z^{(1)}_{t+2}$$z^{(1)}_{t+3}$$y_{t}$$y_{t+1}$$y_{t+2}$$y_{t+3}$$y_{1:t-1}$$y_{1:t-1}$$y_{1:t+1}$$y_{1:t-1}$$y_{1:t}$$y_{1:t+1}$$y_{1:t+2}$

(b) Autoregressive form

Figure 2: Two versions of a hierarchical generative model for sequential data,
$y$, on state space form (top) and autoregressive form (bottom), with latent
variables, $z$. $z^{(2)}$ updates only every other timestep and $z^{(3)}$
updates only every forth.

##  2 Background

In this section we will introduce and give some background on the components
needed to realize the proposed model.

###  2.1 Temporal convolutional network

The temporal Convolutional Network (TCN) is a type of non-linear
autoregressive (NARX) model that makes use of a deep autoregressive neural
network model. Bai et al. (2018) observe that these autoregressive models are
easier to train, less sensitive to hyperparameter selection and in many
circumstances achieves superior results, compared to deterministic state space
models, i.e. LSTMs and RNNs. This can partly be explained by the efficient
implementation, both computationally and parameterwise, of the autoregressive
neural networks here realized with convolutions.

The depth of the TCN is achieved by stacking multiple affine transformations
intertwined with activation functions. Each affine transformation/activation
function pair can be expressed with

$\displaystyle h_{t}=\phi(x_{t-k:t-1}),$ (1)

where $\phi$ combines an affine transformation and a nonlinear activation
function of the regression vector, $x_{t-k:t-1}$. The notation $x_{i:j}$
denotes the sequence of $x$ from $i$ to $j$. An alternative formulation of
this is to view the entire $h_{1:T}$ as a convolution between $\phi$ and the
input $x_{1:T}$. The width of this convolution filter, $k$, corresponds to the
so called the _receptive field_ of the convolution. We use the notation

$\displaystyle h_{1:T}=f(\text{Conv}(x_{1:T})),$ (2)

where Conv stands for the affine transformation and $f$ is an activation
function.

The dilated convolution (Van Den Oord et al., 2016) is an option to the
convolution that sets some of the parameters in the affine transformation to
zero. With a dilation rate $\mu$ only every $\mu$:th element in the affine
transformation is nonzero. By convention $k$ denotes the number of effective
parameters in this transformation. Thus, this increase the receptive field of
the filter to $k\mu$. By stacking several dilated convolutions with
exponentially increasing dilation rate it is possible to produce a model with
exponentially long memory (exponentially large respective field) with a linear
increase in the number of parameters (Yu and Koltun, 2016). The parameters of
the affine transformation are the parameters of the convolutional layer, while
both $k$ and $\mu$ are hyperparameters.

###  2.2 Wavenet

One of the most prominent TCNs is Wavenet (Van Den Oord et al., 2016). Wavenet
is built up of stacked so-called Wavenet blocks, where each block consists of
a residual connection and an identity connection (inspired by ResNet (He et
al., 2016)). Each residual connection in turn consists of a sequence of: a
dilated convolution, a gate-like activation function, and a final $1{\times}1$
convolution. An $1{\times}1$ convolution is a type of convolution that only
operates locally (the filter has size 1). The Wavenet block can be summarized
as

$\displaystyle\text{WavenetBlock}(x_{1:T})$ $\displaystyle=x_{1:T}+$
$\displaystyle\text{Conv}_{1{\times}1}(\tanh($
$\displaystyle\text{Conv}(x_{1:T}))\odot\sigma(\text{Conv}(x_{1:T}))),$ (3)

where $\odot$ denotes element-wise multiplication and $\sigma$ denotes the
sigmoid activation function.

When merging two data sequences with the Wavenet model Van Den Oord et al.
(2016) proposed to consider one of the sequences as the main sequence
($x_{1:T}$) and the other as a conditional sequence ($c_{1:T}$). The two
sequences are concatenated (in feature space) forming a new stream and this is
used as input to the residual network. The identity connection only passes
$x_{1:T}$ forward. Thus, a conditional Wavenet block can be written as,

$\displaystyle\text{CondWavenetBlock}(x_{1:T},c_{1:T})=x_{1:T}+$
$\displaystyle\text{Conv}_{1{\times}1}(\tanh(\text{Conv}([x_{1:T},c_{1:T}]))\odot\sigma(\text{Conv}([x_{1:T},c_{1:T}]))),$
(4)

where $[\cdot]$ denotes concatenation in the feature dimension. When stacking
several such blocks together, all blocks share the same conditioning sequence.
Lai et al. (2018) with Stochastic Wavenet and later Aksan and Hilliges (2019)
extended these models to also include latent variables similar to the proposed
model.

###  2.3 Variational autoencoder

A common approach for latent variables in deep learning is to use variational
autoencoders (VAE) (Rezende et al., 2014; Kingma and Welling, 2014). This
section aims to introduce the notation we use in this paper and for a more in
depth and pedagogical description of the VAE see Kingma and Welling (2019).
The idea of a VAE is to model the distribution of $y$ with a generative model
as

$p(y)=\operatorname*{\mathbb{E}}\big{[}p(y\,|\,z)p(z)\big{]},$ (5)

with

$\displaystyle p(y\,|\,z)$
$\displaystyle=\mathcal{N}(y\,|\,\mu_{y}(z),\Sigma_{y}(z)),$ $\displaystyle
p(z)$ $\displaystyle=\mathcal{N}(z\,|\,0,I),$

where $\mu_{y}(z)$ and $\Sigma_{y}(z)$ are modeled as neural networks.

In the VAE setting the expectation in (5) is estimated with Monte Carlo
samples and (amortized) variational inference. The variational inference is
done by introducing a parameterized approximate posterior distribution, $q$,

$\displaystyle q(z)=\mathcal{N}(z\,|\,\mu_{q}(y),\Sigma_{q}(y)),$

where $\mu_{q}(y)$ and $\Sigma_{q}(y)$ are also modeled as neural networks. It
is also possible to condition the VAE on some arbitrary data, $c$, in which
case the prior is similarly expressed with neural networks,

$\displaystyle p(z\,|\,c)$
$\displaystyle=\mathcal{N}(z\,|\,\mu_{p}(c),\Sigma_{p}(c)).$ (6)

The approximate posterior $q(z)$ and the generative distribution will then
also depend on $c$.

Using this approximate posterior and Jensen’s inequality, we can bound the
log-likelihood from below. This gives us the evidence lower bound (ELBO)

$\mathcal{L}_{\text{ELBO}}(y;\theta)=\mathbb{E}_{q}\big{[}\log
p(y\,|\,z)\big{]}-D_{\text{KL}}\left(q(z)||p(z)\right),$ (7)

where $D_{\text{KL}}$ denotes the Kullback-Lieber divergence and $\theta$ are
the parameters of both the approximative posterior and the generative model.
The expectation here is evaluated with respect to $q$ and it is most commonly
approximated with a single Monte Carlo sample. The ELBO objective is used to
update both the parameters of the generative model and the approximate
posterior.

The variational RNN (VRNN) (Chung et al., 2015), stochastic RNN (SRNN)
(Fraccaro et al., 2016) and STORN (Bayer and Osendorfer, 2014) are all
examples of VAEs that have been adapted to sequential data. The core of these
models is an RNN which is extended with a latent variable. These models
require a backward flow for inference which is implemented by an RNN running
backwards in time.

###  2.4 Hierarchical variational autoencoders

A hierarchical variational autoencoder extends the VAE by introducing
hierarchically stacked latent variables, $\\{z^{(l)}\\}_{l=1}^{L}$, in $L$
layers. The use of stacked latent variables was proposed along with the
original variational autoencoder, though with small improvements compared to
the single layer VAE. The hierarchical models were later improved by Sønderby
et al. (2016) with the so called ladder VAE, where the approximate posterior
is linked to the generative distribution (see Figure 3(a)). This idea has
since then proved to be fertile ground for numerous other architectures, e.g.
ResNet VAE (Kingma et al., 2016), BIVA (Maaløe et al., 2019) and NVAE (Vahdat
and Kautz, 2020).

The approximate posterior’s dependence on $y$ is implemented by extracting a
bottom-up hierarchy of features, $d^{(l)}$, where each successive layer of
features depends on the previous layer, $d^{(l-1)}$. The link between the
generative and the approximate posterior distribution has proven to be
essential for efficient training of these models. Although, the implementation
details for this link varies it can represented as Figure 3(a) in that the
posterior also depends on the extracted features $h^{(l)}$ from the generative
model. The parameters of the approximate posterior are optimized jointly with
the parameters of the generative model.

$d^{(1)}$pt$h^{(1)}$$z^{(1)}$$z^{(1)}$ptKLpt$d^{(2)}$pt$h^{(2)}$$z^{(2)}$$z^{(2)}$ptKLpt$d^{(3)}$pt$h^{(3)}$$z^{(3)}$$z^{(3)}$ptKLpt$y$$y$

(a) Hierarchical VAE. The vertical dashed line was not present in the original
ladder VAE (Sønderby et al., 2016) but it was used in the ResNet VAE (Kingma
et al., 2016).

$d^{(1)}$ptWavenet$h^{(1)}$Delay$z^{(1)}$$z^{(1)}$ptKLpt1x1
CondWavenet$d^{(2)}$ptWavenetDownsample$h^{(2)}$Delay$z^{(2)}$$z^{(2)}$ptKLpt1x1
CondWavenetUpsample1x1
Wavenet$d^{(3)}$ptWavenetDownsample$h^{(3)}$Delay$z^{(3)}$$z^{(3)}$ptKLpt1x1
CondWavenetUpsample1x1 Wavenet$y$$y$

(b) The proposed model, the horizontal dashed line enables the useful
autoregressive relationship.

Figure 3: Approximative posterior model, generative model and common features
for hierarchical VAEs (3(a)) and our model (3(b)). The dotted line encodes the
KL divergence between the prior distribution and the approximate posterior.
Note that during likelihood evaluation $z^{(l)}$ is sampled from the
posterior.

##  3 Model

In our model we combine temporal convolutional networks and striding as
building blocks to instantiate Figure 2(b). We combine this with a
hierarchical VAE where we condition the latent states on the previously
observed data in an autoregressive fashion, as visualized in Figure 4(a).
Analogously to the hierarchical VAE, the model can roughly be split into to a
bottom-up (from data to more abstract features) and a top-down (from abstract
features to predictions) network (Maaløe et al., 2019). The autoregressive
features will also be used as features for the approximate posterior, inspired
by the work of Aksan and Hilliges (2019) which this work extends.

The bottom-up network produces a set of features that is divided into
hierarchical layers. This network also downsamples the signal in between every
layer to capture the multiscale property of the model. To avoid excessive
cluttering we use the notation $d^{(l)}=d^{(l)}_{1:S^{(l)}:T}$ to denote the
bottom-up features which have been downsampled by a total factor of $S^{(l)}$.
The downsample is implemented inside the first WavenetBlock function where the
residual connection simply uses strided convolutions and the identity
connection use the average. We use superscript $D$ to denote Wavenet with
downsampling.

$\displaystyle d^{(1)}$ $\displaystyle=\text{Wavenet}(y),$ (8a) $\displaystyle
d^{(l)}$ $\displaystyle=\text{Wavenet}^{\text{D}}(d^{(l-1)}),\qquad l\neq 1.$
(8b)

The features of the top-down network, $h^{(l)}~{}=~{}h^{(l)}_{1:S^{(l)}:T}$,
follows the hierarchical VAE with an additional conditioning and upsampling.
The conditioning is performed by concatenating the output from the previous
hierarchical layer, which we denote $g^{(l)}=g^{(l)}_{1:S^{(l)}:T}$, with the
delayed features of the bottom-up network, followed by a $1{\times}1$
convolution

$\displaystyle h^{(l)}$
$\displaystyle=\text{Wavenet}_{1{\times}1}([g^{(l)},\text{Delay}(d^{(l)})]),\qquad
l\neq L,$ (9a) $\displaystyle h^{(L)}$
$\displaystyle=\text{Wavenet}_{1{\times}1}(\text{Delay}(d^{(L)})),$ (9b)

where $[\cdot]$ denotes concatenation in feature domain. The Delay function
delays $d^{(L)}$ one step at the current downsample level and prepends zeros,
i.e.

$\displaystyle\text{Delay}(d_{1:S:T})=[0,d_{1:S:T-S}],$ (10)

here the concatenation is done in time domain. Finally
$\text{Wavenet}_{1{\times}1}$ denotes a Wavenet where all the dilated
convolutions are replaced with $1{\times}1$ convolutions.

The prior distribution for the latent variable is given by

$\begin{split}p(z^{(l)}\,|\,h^{(l)})=&\\\
\mathcal{N}\bigg{(}z^{(l)}&\bigg{|}\begin{split}&\mu=\text{Conv}_{1{\times}1}(h^{(l)}),\\\
&\Sigma=\text{diag}(\text{softplus}(\text{Conv}_{1{\times}1}(h^{(l)}))^{2})\end{split}\bigg{)},\end{split}$
(11)

Finally, the output of the top-down layer is calculated as,

$\displaystyle g^{(0)}$
$\displaystyle=\text{CondWavenet}_{1{\times}1}(h^{(1)},z^{(1)}),$ (12a)
$\displaystyle g^{(l-1)}$
$\displaystyle=\text{CondWavenet}^{\text{U}}_{1{\times}1}(h^{(l)},z^{(l)}),\,\,\,l\neq
1,$ (12b)

thus, the latent variables are used as a conditioning sequence in the Wavenet.
The Upsampling here, $\text{CondWavenet}^{\text{U}}$, is similar to the
downsample made as a part of the Wavenet. The identity connection is hear
upsampled with nearest neighbor and the residual with strided transpose
convolutions. The predictive distribution is chosen as
$p(y\,|\,\text{Conv}_{1{\times}1}(g^{(0)}))$. The shape of the predictive
distribution is chosen with respect to the problem – Gaussian for real valued
data and Bernoulli for binary data or mixtures thereof.

The approximate posterior is expressed using the previously defined features
$d^{(l)}$ and $h^{(l)}$ as

$\begin{split}q(z^{(l)})=&\\\
\mathcal{N}\bigg{(}&z^{(l)}\bigg{|}\begin{split}&\mu=\text{Conv}_{1{\times}1}([h^{(l)},d^{(l)}]),\\\
&\Sigma=\text{diag}(\text{softplus}(\text{Conv}_{1{\times}1}([h^{(l)},d^{(l)}])^{2})\end{split}\bigg{)}.\end{split}$
(13)

We do not calculate $q$ as a combination of $p$ and another normal
distribution, as advocated by Sønderby et al. (2016) and Aksan and Hilliges
(2019). Instead $q$ is parameterized independently of the common features,
similar to Kingma et al. (2016).

Figures 4(a) and 4(b) visualize the model in an alternative fashion compared
to Figure 3(b). Figure 4(a) resolves the temporal dimension and Figure 4(b)
visualizes a hierarchical layer locally in time. The parameters of this model
(both the generative model and the approximate posterior) are the parameters
for the Wavenets, the Convolutional layers, the Upsample layers and the
Downsample layers.

Although the stochastic TCN (STCN) (Aksan and Hilliges, 2019) is very related
to our model, there are some distinctions. We generalize the STCN with the
hierarchical VAE in favor of the Ladder VAE and also extend it by
incorporating multiple timescales. Aksan and Hilliges (2019) argue that their
model expresses multiple timescales but we disagree since it is missing the
downsampling and upsampling steps. Indeed, the deeper features of the STCN
gathers information from a larger receptive field. However, these features do
only affect the current time step and are not shared between different
timesteps.

pt$z^{(3)}_{t}$$z^{(2)}_{t}$$z^{(2)}_{t+2}$$z^{(1)}_{t}$$z^{(1)}_{t+1}$$z^{(1)}_{t+2}$$z^{(1)}_{t+3}$$y_{t}$$y_{t+1}$$y_{t+2}$$y_{t+3}$$d^{(3)}_{t-4}$$d^{(2)}_{t-2}$$d^{(2))}_{t}$$d^{(1)}_{t-1}$$d^{(1)}_{t}$$d^{(1)}_{t+1}$$d^{(1)}_{t+2}$

(a) The hierarchical structure of the generative model for a sample with four
timesteps. This structure repeats to the right for longer sequences. The gray
arrows corresponds to the bottom-up paths. Note that the $d$s are delayed with
$S^{(l)}$.

$\cdots\text{Downsample}\cdots$$d_{t-2k}^{(l-1)}$$d_{t+1}^{(l-1)}$$d^{(l)}_{t}$Wavenetpt$h^{(l)}_{t+2}$$z^{(l)}_{t+2}$ptptUpsampleptptpt$g_{t+2}^{(l-1)}$$g_{t+3}^{(l-1)}$$g_{t+2}^{(l)}$$d_{t}^{(l)}$

(b) A close-up of the dashed square in Figure 4(a) for the generative model.
Here $r$ is the receptive field of the Wavenet model. In this instance
$S^{(l)}=2$.

Figure 4: An alternative description of the generative part of the model in
Figure 3(b) where the temporal dimension is explicitly shown.

##  4 Experiments

We conduct experiments on two hierarchically structured sequential datasets,
raw waveform speech (Blizzard) and handwritten text (IAM-OnDB). For all
experiments we use linear free bits (see Appendix A) and the Adam optimizer
with an exponentially decaying learning rate of $5\cdot 10^{-4}$. For free
bits we used a scheme to smoothly decay the threshold of free bits
exponentially over the course of training at a rate of halving the threshold
every $300\thinspace 000$ iterations. We train for $100$ epochs and do not see
any tendencies of overfitting on the training data. Finally, each bottom-up
Wavenet consists of four Wavenet blocks and each top-down constsits of two
Wavenet blocks. The number of Wavenet blocks are thus the same in the bottom-
up and the top-down networks. For exact hyperparameter settings and additional
implementation details we refer to the appendix and
code111https://github.com/carl-andersson/MS-STCN.

Blizzard is a raw waveform speech dataset that originates from the 2013
Blizzard challenge (Karaiskos and Vasilis, 2013). The raw signal is split into
$0.5$ second ($8\thinspace 000$ samples) long sequences. We experiment by
dividing these sequences into frames of $200$ samples down to $2$ to really
enable the long memory requirements and hierarchical decomposition of the
data. Previous research on this dataset used frames of $200$ samples.
Additionally, we remove some artifacts from the mp3 decoding that appear
before and after every sequence. The remaining preprocessing is performed
according to Aksan and Hilliges (2019).

The smaller framesize increases the computational complexity of the model
compared to that of the STCN if we were to use the same number of filters. We
compensate for this by decreasing the number of filters in our model for
layers with small $l$ compared to those of the STCN. This way the
computational complexity is approximately the same for the all of proposed
models as for the (moderate) STCN, while at the same time using significantly
fewer parameters. In the case of frame width of 200 we keep all the filters
the same size which to be comparable with STCN.

The ELBO evaluated on test data for the Blizzard dataset is presented in Table
1. Our model without striding is essentially a reimplementation of STCN using
the model setup described in this paper. As a control experiment we also
report the result with our model without the stochastic latent variables (i.e.
only deterministic autoregressive). However, this control showed instabilities
and eventually diverged so we report the best fit we got before it diverged.
This divergence might be an indication that the latent variables are somehow
regularizing the model. In Table 2 we also report the amount of information
that is used in the KL-term for the different layers. Results show that KL
units on all levels are being used.

IAM-OnDB (Liwicki and Bunke, 2005) is a handwritten text dataset where every
sample consists of the pen position ($\mathbb{R}^{2}$) and a binary variable
indicating whether the pen was lifted after this position or not. We use the
same preprocessing as Chung et al. (2015) and Aksan and Hilliges (2019).
However, we noted that the discretization of the pen position in some cases
led to models that greatly overestimated the likelihood fit to the data. To
cope with these discretization artifacts we added artificial uniform
measurement noise to the pen position.

Samples from the model trained on IAM-OmDB (Figure 1) shows complex sequences
were it is even possible to recognize a few words. However, the added noise
means that the likelihood estimates of the model are not easily comparable
with previous results on this dataset.

Table 1: Average ELBO per sequence comparison for speech with varying frame width (FW) and with and without down-/upsampling (multiscale).  Model | Blizzard
---|---
RNN (GMM) (Chung et al., 2015) | $7413$
Wavenet (GMM) Aksan and Hilliges (2019) | $8190$
VRNN (GMM) (Chung et al., 2015) | $\approx 9392$
SRNN (Fraccaro et al., 2016) | $\geq 11\,991$
Z-Forcing (Goyal et al., 2017) | $\geq 15\,430$
STCN (Moderate) (Aksan and Hilliges, 2019) | $\geq 16\,288$
STCN (Large) (Aksan and Hilliges, 2019) | $\geq 17\,670$
Our model, FW=200 | $\geq 13\,645$
Our model, FW=200 w/o multiscale | $\geq 15\,746$
Our model, FW=25 | $\geq 17\,060$
Our model, FW=25 w/o multiscale | $\geq 16\,562$
Our model, FW=2 | $\geq 18\,809$
Our model, FW=2 w/o multiscale | $\geq 18\,372$
Our model, FW=2 w/o latent variables | $18\,523$

##  5 Conclusion and future work

In this paper we argue that a deep learning model is not that different from a
structured autoregressive model. However, the philosophy behind the model and
the problems it is applied to, are significantly different from a typical
system identification model. The model we proposed has shown to model
hierarchically structured data very well while at the same time being more
parameter efficient compared to existing models. We hypothesize that the
proposed model exploit the hierarchical nature of the sequential data similar
to how the CNN exploits the hierarchical structure in image data.

The proposed model is in principle not limited to one-dimensional temporal
data. Similar to PixelVAE (Gulrajani et al., 2017), our model can be applied
to images or even video. Investigating whether the promising results of this
model transfers into other regimes of hierarchically structured data is an
interesting avenue for future work. Another area for further experiments is to
apply the model to other hierarchically structured data discussed in the
introduction, one such example is pose and orientation estimation from
inertial sensors (Kok et al., 2017).

Table 2: Average KL-divergence per layer on the Blizzard dataset where KL1 corresponds to the layer closest from the data.  Model | KL1 | KL2 | KL3 | KL4 | KL5
---|---|---|---|---|---
FW=200 | $0.584\,920\,180\,000\,505\,2$ | $7.990\,862\,283\,759\,18$ | $4.451\,699\,510\,602\,553$ | $3.931\,846\,768\,133\,616$ | $74.692\,176\,265\,760\,35$
FW=200 w/o multiscale | $48.639\,276\,380\,958\,35$ | $33.656\,451\,782\,743\,52$ | $47.963\,004\,908\,853\,83$ | $146.007\,384\,007\,717\,18$ | $142.979\,792\,928\,635\,62$
FW=25 | $235.4$ | $45.289\,409\,35$ | $44.335\,122\,93$ | $66.625\,844\,68$ | $101.700\,344\,07$
FW=25 w/o multiscale | $320.934\,691\,483\,796\,03$ | $126.599\,087\,071\,288\,65$ | $488.131\,012\,366\,204\,63$ | $542.951\,341\,801\,324\,1$ | $180.425\,208\,261\,803\,87$
FW=2 | $629.181\,338\,474\,073\,9$ | $48.945\,017\,299\,496\,2$ | $45.187\,340\,858\,576\,84$ | $35.259\,783\,210\,851\,42$ | $57.979\,644\,418\,383\,87$
FW=2 w/o multiscale | $203.507\,759\,389\,451\,8$ | $210.291\,654\,083\,148\,7$ | $158.277\,840\,564\,747\,66$ | $209.417\,462\,246\,597\,6$ | $50.336\,778\,243\,226\,17$

## Acknowledgment

We thank Joakim Jaldén’s group at KTH, Royal Institute of Technology,
Stockholm, for lending computational power.

## References

  * Aksan and Hilliges (2019) Aksan, E. and Hilliges, O. (2019).  STCN: Stochastic Temporal Convolutional Networks.  In _International Conference on Learning Representations_.
  * Bai et al. (2018) Bai, S., Kolter, J.Z., and Koltun, V. (2018).  An Empirical Evaluation of Generic Convolutional and Recurrent Networks for Sequence Modeling.  URL http://arxiv.org/abs/1803.01271.
  * Bayer and Osendorfer (2014) Bayer, J. and Osendorfer, C. (2014).  Learning Stochastic Recurrent Networks.  URL http://arxiv.org/abs/1411.7610.
  * Bowman et al. (2016) Bowman, S.R., Vilnis, L., Vinyals, O., Dai, A., Jozefowicz, R., and Bengio, S. (2016).  Generating Sentences from a Continuous Space.  In _Conference on Computational Natural Language Learning_.
  * Chung et al. (2016) Chung, J., Ahn, S., and Bengio, Y. (2016).  Hierarchical Multiscale Recurrent Neural Networks.  URL http://arxiv.org/abs/1609.01704.
  * Chung et al. (2015) Chung, J., Kastner, K., Dinh, L., Goel, K., Courville, A., and Bengio, Y. (2015).  A Recurrent Latent Variable Model for Sequential Data.  In _Advances in Neural Information Processing Systems_.
  * Fraccaro et al. (2016) Fraccaro, M., Sønderby, S.K., Paquet, U., and Winther, O. (2016).  Sequential Neural Models with Stochastic Layers.  In _Advances in Neural Information Processing Systems_. Curran Associates, Inc.
  * Goyal et al. (2017) Goyal, A., Sordoni, A., Maluuba, M., Côté, M.A., Rosemary, N., Mila, K., Montréal, P., and Bengio, Y. (2017).  Z-Forcing: Training Stochastic Recurrent Networks.  URL https://arxiv.org/abs/1711.05411.
  * Gulrajani et al. (2017) Gulrajani, I., Kumar, K., Ahmed, F., Taiga, A.A., Visin, F., Vazquez, D., and Courville, A. (2017).  PixelVAE: A Latent Variable Model for Natural Images.  In _International Conference on Learning Representations_.
  * He et al. (2016) He, K., Zhang, X., Ren, S., and Sun, J. (2016).  Deep Residual Learning for Image Recognition.  In _Conference on Computer Vision and Pattern Recognition_.
  * Hochreiter and Urgen Schmidhuber (1997) Hochreiter, S. and Urgen Schmidhuber, J. (1997).  Long Short-Term Memory.  _Neural Computation_ , 9.
  * Karaiskos and Vasilis (2013) Karaiskos, S.K. and Vasilis (2013).  The Blizzard Challenge 2013.  In _Blizzard Challenge workshop_.
  * Kingma and Welling (2014) Kingma, D.P. and Welling, M. (2014).  Auto-Encoding Variational Bayes.  In _International Conference on Learning Representations_.
  * Kingma and Welling (2019) Kingma, D.P. and Welling, M. (2019).  An Introduction to Variational Autoencoders.  _Foundations and Trends in Machine Learning_.
  * Kingma et al. (2016) Kingma, D.P., Salimans, T., Jozefowicz, R., Chen, X., Sutskever, I., and Welling, M. (2016).  Improved Variational Inference with Inverse Autoregressive Flow.  In _Advances in Neural Information Processing Systems_.
  * Kok et al. (2017) Kok, M., Hol, J.D., and Schön, T.B. (2017).  Using inertial sensors for position and orientation estimation.  _Foundations and Trends® in Signal Processing_.
  * Koutník et al. (2014) Koutník, J., Greff, K., Gomez, F., and Schmidhuber, J. (2014).  A clockwork RNN.  In _International Conference on Machine Learning_.
  * Krizhevsky et al. (2012) Krizhevsky, A., Sutskever, I., and Hinton, G.E. (2012).  ImageNet Classification with Deep Convolutional Neural Networks.  In _Advances in Neural Information Processing Systems_.
  * Lai et al. (2018) Lai, G., Li, B., Zheng, G., and Yang, Y. (2018).  Stochastic WaveNet: A Generative Latent Variable Model for Sequential Data.  URL https://arxiv.org/abs/1806.06116.
  * LeCun et al. (2015) LeCun, Y., Bengio, Y., Hinton, G., Y., L., Y., B., and G., H. (2015).  Deep learning.  _Nature_ , 521.
  * Liwicki and Bunke (2005) Liwicki, M. and Bunke, H. (2005).  IAM-OnDB - an On-Line English Sentence Database Acquired from Handwritten Text on a Whiteboard.  In _International Conference on Document Analysis and Recognition_.
  * Maaløe et al. (2019) Maaløe, L., Fraccaro, M., Liévin, V., and Winther, O. (2019).  BIVA: A Very Deep Hierarchy of Latent Variables for Generative Modeling.  URL https://arxiv.org/abs/1902.02102.
  * Paszke et al. (2017) Paszke, A., Gross, S., Chintala, S., Chanan, G., Yang, E., DeVito, Z., Lin, Z., Desmaison, A., Antiga, L., and Lerer, A. (2017).  Automatic differentiation in pytorch.  In _Advances in Neural Information Processing Systems, Workshop Autodiff_.
  * Rezende et al. (2014) Rezende, D.J., Mohamed, S., and Wierstra, D. (2014).  Stochastic Backpropagation and Approximate Inference in Deep Generative Models.  In _International Conference on Machine Learning_.
  * Schoukens (2020) Schoukens, M. (2020).  Nonlinear benchmark website.  URL http://www.nonlinearbenchmark.org.
  * Sjöberg et al. (1995) Sjöberg, J., Zhang, Q., Ljung, L., Benveniste, A., Delyon, B., Glorennec, P.Y., Hjalmarsson, H., and Juditsky, A. (1995).  Nonlinear black-box modeling in system identification: a unified overview.  _Automatica_ , 31.
  * Sønderby et al. (2016) Sønderby, C.K., Raiko, T., Maaløe, L., Sønderby, S.K., and Winther, O. (2016).  Ladder Variational Autoencoders.  In _Advances in Neural Information Processing Systems_.
  * Vahdat and Kautz (2020) Vahdat, A. and Kautz, J. (2020).  NVAE: A deep hierarchical variational autoencoder.  In _Neural Information Processing Systems_.
  * Van Den Oord et al. (2016) Van Den Oord, A., Dieleman, S., Zen, H., Simonyan, K., Vinyals, O., Graves, A., Kalchbrenner, N., Senior, A., and Kavukcuoglu, K. (2016).  WaveNet: A Generative Model for Raw Audio.  URL https://arxiv.org/abs/1609.03499.
  * Yu and Koltun (2016) Yu, F. and Koltun, V. (2016).  Multi-scale context aggregation by dilated convolutions.  In _International Conference on Learning Representations_.

##  Appendix A Hyperparameters and Implementation details

Code is implemented in Pytorch (Paszke et al., 2017). All experiments are done
with a Titan Xp except the final experiment on the Blizzard dataset which was
done with a Tesla V100.

###  A.1 Hyperparmeters

This section specifies the hyperparameters for the different models used. All
models consists of 5 hierarchical layers. However the number of filters are
chosen to match that the computational load/number of parameters in STCN.

Blizzard: The Wavenets in the bottom-up network consists of 5 Wavenet blocks
with kernel size 2 and the Wavenets in the top-down network consts of 2
Wavenet Blocks each. The dilation rates for the bottom-up networks where
chosen as $[1,2,4,8,16]$ when not using downsampling and $[]1,1,1,1,1]$
otherwise. An exception to this is for frame widths 25 and 2 with downsampling
as the topmost bottom-up layer then uses $[1,2,4,8,16]$

Frame width 200 and multiscale 256 filters in all convolutions with 64 latent
variables for each VAE layer and timestep. The downsampling is $[1,2,2,2,1]$
for the VAE layers respectively and the predictive distribution is chosen as a
mixture of 10 Gaussians.

Frame width 200 and NO multiscale 256 filters in all convolutions with 64
latent variables for each VAE layer and timestep. The downsampling rate is
$[1,1,1,1,1]$ for the VAE layers respectively and the predictive distribution
is chosen as a mixture of 10 Gaussians.

Frame width 25 and multiscale $[96,128,192,256,256]$ filters the convolutions
in respective VAE layer and $[24,32,48,64,64]$ latent variables for each VAE
layer and timestep. The downsampling rate is $[1,2,2,2,1]$ for the VAE layers
respectively and the predictive distribution is chosen as a mixture of 10
Gaussians.

Frame width 25 and NO multiscale 96 filters for all the convolutions in
respective VAE layer and 24 latent variables for each VAE layer and timestep.
The downsampling rate is $[1,1,1,1,1]$ for the VAE layers respectively and the
predictive distribution is chosen as a mixture of 10 Gaussians.

Frame width 2 and multiscale $[26,58,128,256,256]$ filters the convolutions in
respective VAE layer and $[1,5,25,64,64]$ latent variables for each VAE layer
and timestep. The downsampling rate is $[1,5,5,4,1]$ for the VAE layers
respectively and the predictive distribution is chosen as a mixture of 2
Gaussians.

Frame width 2 and NO multiscale 26 filters for all the convolutions in
respective VAE layer and 1 latent variables for each VAE layer and timestep.
The downsampling rate is $[1,1,1,1,1]$ for the VAE layers respectively and the
predictive distribution is chosen as a mixture of 2 Gaussians.

The number of filter for each examples is chosen to keep the computational
complexity equal to that of the moderate STCN.

IAM-OnDB The Wavenets in the bottom-up network consists of 4 Wavenet blocks
with kernel size 2 and the Wavenets in the top-down network consts of 2
Wavenet Blocks each. The dilation rates for the bottom-up networks where
chosen as $[1,1,1,1,1]$ except for the topmost layer where we used
$[1,2,4,8]$. 64 filters in all convolutions with $[8,8,8,4,2]$ latent
variables for VAE layer respectively and timestep. The downsampling is
$[1,2,2,2,1]$ for the VAE layers respectively and the predictive distribution
is chosen as a mixture of 10 Gaussians and a Bernoulli for the pen up event.

###  A.2 Predictive distributions

The final distribution is parameterized as
$p(y\,|\,z^{(1)})=p(y\,|\,\text{Conv}_{1{\times}1}(g^{(0)}))$. Let us denote
the features $\text{Conv}_{1{\times}1}(g^{(0)})$ with $\eta$ for some clarity,
we also denote the dimension of $y$ with $p$.

#### Gaussian

We parameterize the final prediction the same way as the latent variables. The
features $\eta$ is split up into mean features, $\eta_{\mu}$ and standard
deviation features, $\eta_{\sigma}$. In addition to what is mentioned in the
paper we clamp the standard deviation of the Gaussian distribution between
$10^{-3}$ and $5$, both for the final prediction and the latent variables. The
distribution is thus given of

$\displaystyle
p(y\,|\,z^{(1)})=\mathcal{N}(y\,|\,\eta_{\mu},\text{softmax}(\text{Clamp}(\eta_{\sigma},10^{-3},5.0))^{2}),$
(14)

and the dimension of $\eta$ must thus be $2p$.

#### Bernoulli

For the Bernoulli distribution we simply parameterize with the logits which is
equal to $\eta$. The dimension for $\eta$ is thus $p$.

#### Mixtures

For mixture distributions with $n_{c}$ components we split $\eta$ into $n_{c}$
weight parameters and the rest to parameterize the $n_{c}$ base distributions.
The weight parameters are interpreted as the logit weights. The dimension for
$\eta$ in this case is $2pn_{c}+n_{c}$ for the Gaussian mixture model and
$pn_{c}+n_{c}$ for the Bernoulli mixture model.

###  A.3 Posterior collapse

Without any additional tricks, the hierarchical VAE suffers from _posterior
collapse_ , a state during the optimization where the model gets stuck in a
local minimum and the KL units do not transfer any information (Sønderby et
al., 2016; Kingma et al., 2016; Goyal et al., 2017; Bowman et al., 2016). This
state can be circumvented through a number of different strategies, e.g. _KL-
annealing_ (Sønderby et al., 2016; Bowman et al., 2016), _free bits_ (Kingma
et al., 2016) or auxiliary loss functions as in Z-forcing (Goyal et al.,
2017).

To avoid posterior collapse we use an alternative to free bits (Kingma et al.,
2016). The ordinary free bits sets a minimum value of how low the KL units can
be minimized by altering the ELBO objective Equation 7 to

$\mathcal{L}_{\text{FB}}=\operatorname*{\mathbb{E}}_{q}\log
p(y\,|\,z)-\sum_{i}\text{maximum}\left(D_{\text{KL}}\left(q(z_{i})||p(z_{i})\right),\lambda\right),$
(15)

where $\lambda$ denotes the free bits threshold and $i$ iterates through all
latent variables, $z$. To smoothen the transition from where the threshold is
active to where it is not, we propose a linear smoothing of the free bits
objective. Hence, we define linear free bits

$\begin{split}\mathcal{L}&{}_{\text{LFB}}=\operatorname*{\mathbb{E}}_{q}\log
p(y\,|\,z)-\\\
&\sum_{i}\Psi\left(\text{min}\left(\frac{D_{\text{KL}}\left(q(z_{i})||p(z_{i})\right)}{\lambda},1\right)\right)D_{\text{KL}}\left(q(z_{i})||p(z_{i})\right),\end{split}$
(16)

where $\Psi(\cdot)$ is a stop gradient function that ensures that no gradient
information is passed through the function in the backward pass while acting
as the identity function in the forward pass.


[a,b]Rui Zhang

# Strange and charm contributions to nucleon charges and moments111 The work
of RZ is supported by the US National Science Foundation under grant PHY
1653405 “CAREER: Constraining Parton Distribution Functions for New-Physics
Searches”. We thank MILC Collaboration for sharing the lattices used to
perform this study. The LQCD calculations were performed using the Chroma
software suite [1] with the multigrid solver algorithm [2, 3]. This research
used resources of the National Energy Research Scientific Computing Center, a
DOE Office of Science User Facility supported by the Office of Science of the
U.S. Department of Energy under Contract No. DE-AC02-05CH11231 through ERCAP;
the Extreme Science and Engineering Discovery Environment (XSEDE), which is
supported by National Science Foundation grant number ACI-1548562; and
facilities of the USQCD Collaboration, which are funded by the Office of
Science of the U.S. Department of Energy, Extreme Science and Engineering
Discovery Environment (XSEDE), which is supported by National Science
Foundation grant number ACI-1548562.

Tanmoy Bhattacharya  Rajan Gupta  Huey-Wen Lin  Santanu Mondal  Sungwoo Park
Boram Yoon

###### Abstract

We present preliminary results for strange and charm contributions to nucleon
charges and moments. The scalar, axial and tensor charges, and unpolarized
first moments are calculated using clover-on-HISQ formulation and cover four
lattice spacings, $a=\\{0.06,0.09,0.12,0.15\\}$ fm, and three pion masses,
$M_{\pi}=\\{310,220,130\\}$ MeV. The renormalization factors are calculated
nonperturbatively using the RI-sMOM scheme. We carry out a chiral and
continuum extrapolation to obtain physical results.

##  1 Introduction

Nucleon charges and moments are important quantities to study for the
elucidation of the nucleon structure. The scalar charge plays an important
role in dark matter search, and gives the pion-nucleon sigma term. The axial
charge is the contribution of quark spins to the proton spin, and enters in
the spin-dependent dark matter cross-section. The tensor charge gives the
contribution of the quark Electric Dipole Moment (EDM) to the neutron EDM. The
first unpolarized moment gives the momentum fraction of the quark, an
important quantity to describe the quark PDF especially for sea quarks. These
quantities for the isovector ($u-d$) case have been calculated with high
precision using lattice QCD. As the precision of light quark calculation
improves, the contributions from heavier quark flavors, strange and charm,
also need to be precisely determined.

All operators for charge and moment calculations are quark bilinears,
$O=\bar{q}\Gamma q$. Specifically, for the charge operators,
$\Gamma=\mathbf{1}$ for the scalar charge, $\Gamma=\gamma_{i}\gamma_{5}$ for
the axial charge, $\Gamma=\sigma_{ij}$ for the tensor charge, are local
operators. On the other hand,
$\Gamma=D_{4}\gamma_{4}-\tfrac{1}{3}\sum_{i}D_{i}\gamma_{i}$ is a non-local
operator for the first unpolarized moment, which gives the momentum fraction
carried by the quark.

For strange and charm quarks, the contribution only comes from disconnected
diagrams, ie, the correlation of the self contraction of the quark bilinear
operator $\bar{q}\Gamma q$ (forming a quark loop) with the nucleon, as shown
in Fig. 1. The details of disconnected calculations can be found in Ref. [4].

Figure 1: Left: Disconnected diagram on lattice. Connected NPR diagrams
(Middle) for all operators and mixing diagrams (right) for the scalar charge.

##  2 Lattice Setup, Matrix Elements and Renormalization

We calculated 2- and 3-point functions on six ensembles generated with 2+1+1
flavors of highly improved staggered quarks (HISQ) [5] by the MILC
collaboration [6]. They cover four lattice spacings,
$a\approx\\{0.06,0.09,0.12,0.15\\}$ fm, and 3 pion masses of
$M_{\pi}\approx\\{130,220,310\\}$ MeV. To reduce statistical noise, the
lattices are hypercubic (HYP) smeared [7] before the calculation of the quark
propagators. The details of the ensembles are summarized in Table 1.

a (fm) |  $M_{\pi}$ (MeV) | $L^{3}\times T$ | $N_{\text{conf}}^{c}$ | $N_{\text{conf}}^{s}$ | $N_{\text{src}}^{c}$ | $N_{\text{src}}^{s}$ | $\frac{N_{\text{src}}^{\text{2pt}}}{\text{cfg}}$
---|---|---|---|---|---|---|---
0.06 | 320(2) | $48^{3}\times 144$ | 469 | 87 | 4000 | 4000 | 32
0.09 | 313(3) | $32^{3}\times 96$ | 633 | 889 | 6000 | 6000 | 32
0.09 | 138(1) | $64^{3}\times 96$ | 721 | 310 | 4000 | 4000 | 32
0.12 | 310(3) | $24^{3}\times 64$ | 983 | 897 | 4000 | 8000 | 32
0.12 | 228(2) | $32^{3}\times 64$ | 902 | 869 | 5000 | 5000 | 32
0.15 | 321(4) | $16^{3}\times 48$ | n/a | 1795 | n/a | 2000 | 32
Table 1: Parameters of the 2+1+1-flavor HISQ ensembles used for the strange
and charm contributions to the nucleon charges and moments. $N_{\text{conf}}$
gives the number of gauge configurations, $N_{\text{src}}$ is the number of
the noise source measurements to approximate the disconnected loop
propagators, and $\frac{N_{\text{src}}^{\text{2pt}}}{cfg}$ is the number of
the two-point sources on each configuration.

We computed the nucleon two point correlators of the nucleon interpolating
operator $\chi$ with different projections $\mathcal{P}$:

$C^{\mathcal{P}}_{\text{2pt}}(\tau)=\langle
0|\mathcal{P}\chi(\tau)\bar{\chi}(0)|0\rangle,$ (1)

where $\mathcal{P}^{\text{unp}}=(1+\gamma_{4})/2$ for unpolarized projection,
and $\mathcal{P}^{\text{pol}}=(1+\gamma_{4})(1+i\gamma_{5}\gamma_{i})/2$ for
polarized projection. Then we compute the disconnected loops:

$C^{O}_{\text{loop}}=\sum_{n}\Tr\left[\hat{O}S(n)\right]$ (2)

These two are combined together to obtain the disconnected contribution.

$\displaystyle C^{\text{O}}_{\text{3pt}}(t,\tau)=$
$\displaystyle\left\langle\left(C^{\mathcal{P}\phantom{O}}_{\text{2pt}}(\tau)-\langle
C^{\mathcal{P}\phantom{O}}_{\text{2pt}}(\tau)\rangle\right)\left(C^{O}_{\text{loop}}(t)-\langle
C^{O}_{\text{loop}}(t)\rangle\right)\right\rangle$

We make 2-state fits to the correlators using their spectral decomposition:

$C^{\mathcal{P}^{unp}}_{\text{2pt}}(\tau)=|A_{0}|^{2}e^{-E_{0}\tau}+|A_{1}|^{2}e^{-E_{1}\tau}$
(3) $\displaystyle C^{O}_{\text{3pt}}(t,\tau)=\left(|A_{0}|^{2}\langle
0|O|0\rangle e^{-E_{0}t}e^{-E_{0}(\tau-t)}+A_{0}A_{1}\langle 1|O|0\rangle
e^{-E_{0}t}e^{-E_{1}(\tau-t)}+\right.$ $\displaystyle\left.A_{1}A_{0}\langle
0|O|1\rangle e^{-E_{1}t}e^{-E_{0}(\tau-t)}+|A_{1}|^{2}\langle 1|O|1\rangle
e^{-E_{1}t}e^{-E_{1}(\tau-t)}\right)$ (4)

to extract the ground state matrix elements $\langle 0|O|0\rangle$. The fit
results are displayed in Fig. 2 along with the data for the ratio of 3-point
to 2-point functions,
$R(t,\tau)=C^{O}_{\text{3pt}}(t,\tau)/C^{\mathcal{P}^{unp}}_{\text{2pt}}(\tau)$.
The fit results are consistent with the data points.

Figure 2: Ratio plots from lattice data (colored markers) and the
reconstructed ratio from 2-state fit results (colored bands) on the a12m310
ensemble. The gray bands label the ground state matrix elements.

To renormalize the matrix elements, we used the regularization-independent
symmetric momentum subtraction (RI-sMOM) scheme. In this scheme, we choose
symmetric momentum configurations for the external quark state and the
momentum transfer $p_{1}^{2}=p_{2}^{2}=(p_{2}-p_{1})^{2}=\mu_{R}^{2}$, as
shown in the middle panel of Fig. 1. They are then matched to $\overline{MS}$
scheme using 2-loop perturbation theory, and run to $\mu=2$ GeV. The
$\mu_{R}$-dependence is then removed by fitting the $Z$-factors to
$Z_{\Gamma}(\mu)=Z_{\Gamma}+c_{1}\mu_{R}^{2}+c_{2}\mu_{R}^{4}$. The
renormalization for the scalar charge and the momentum fraction is more
complicated. The momentum fraction for quarks mixes with other flavors as well
as gluons, which is not included in our work so far. The scalar charge also
mixes with other flavors, as shown in the right of Fig. 1, and we have
included the mixing between strange quark and light iso-scalar combination,
but not for the charm quark so far. The two diagrams for calculating $Z_{S}$
are done for all the flavor combinations as the connected ($c_{ff^{\prime}}$)
and disconnected ($d_{ff^{\prime}}$) parts:

$c_{ff^{\prime}}\equiv\frac{\delta_{ff^{\prime}}}{12Z^{f}_{\psi}}\Tr\left[\langle
f|O^{f^{\prime}}|f\rangle_{\text{conn}}\right],\,d_{ff^{\prime}}\equiv\frac{-1}{12Z^{f}_{\psi}}\Tr\left[\langle
f|O^{f^{\prime}}|f\rangle_{\text{disc}}\right].$ (5)

Using these, we renormalize the bare charges in the following way:

$\begin{pmatrix}g_{S}^{R,u+d}\\\ g_{S}^{R,s}\\\
g_{S}^{R,c}\end{pmatrix}=\begin{pmatrix}c_{l}-2d_{ll}&-2d_{sl}&-2d_{cl}\\\
-d_{ls}&c_{s}-d_{ss}&-d_{cs}\\\
-d_{lc}&-d_{sc}&c_{c}-d_{cc}\end{pmatrix}^{-1}\cdot\begin{pmatrix}g_{S}^{\text{bare,u+d}}\\\
g_{S}^{\text{bare,s}}\\\ g_{S}^{\text{bare,c}}\end{pmatrix}$ (6)

Note that $d_{ff^{\prime}}$ is sensitive to the flavor in the quark loop,
$f^{\prime}$, but insensitive to the external quark $f$.

##  3 Chiral and Continuum Extrapolation

After renormalization, we obtain the charges and momentum fraction on
different ensembles. The results for strange quark are shown in Fig. 3, and
those for charm quark are shown in Fig. 4. We extrapolate these quantities to
the physical point by fitting to the ansatz:

$g(M_{\pi},a)=g^{\text{phy}}(1+c_{i}a^{i}+d_{2}(M_{\pi}^{2}-(M_{\pi}^{\text{phy}})^{2})),\qquad
i=1,2,$ (7)

for linear ($i=1$) or quadratic ($i=2$) dependence on the lattice spacing.

Figure 3: Renormalized charges and momentum fraction for the strange quark.

Figure 4: Renormalized charges and momentum fraction for the charm quark.

The extrapolations linear (quadratic) in lattice spacing are shown as magenta
(purple) bands in Fig. 5 and Fig. 6 for strange quark results. We notice that
the $M_{\pi}$ dependence is non-negligible while the $a$ dependence in
$g_{A}$, $g_{T}$, $O_{44}^{S}$ is not clear. The $a$ dependence is large in
$g_{S}$. The linear and quadratic fits have similar $\chi^{2}$ and the two
results are consistent within 1$\sigma$ except for $g_{S}$. We are working on
improving the statistics on finer lattices to better constrain the
extrapolation.

Figure 5: Chiral and continuum extrapolation of strange contribution to
charges and momentum fraction using Eq. (7), and shown as a function of
$M_{\pi}^{2}$ with $a=0$.

Figure 6: Chiral and continuum extrapolation of strange contribution to
charges and momentum fraction using Eq. (7), and shown as a function of
lattice spacing with $M_{\pi}$ set to $135$ MeV.

The extrapolations linear (quadratic) in lattice spacing are shown as magenta
(purple) bands in Fig. 7 and Fig. 8 for charm quark results. Unlike the
strange results, the $M_{\pi}$ dependence is not significant for charm. The
$a$ and $a^{2}$ extrapolations are consistent within 1$\sigma$. We also note
that the scalar charge reduces under continuum extrapolation, while the
momentum fraction becomes larger under continuum extrapolation.

Figure 7: Chiral and continuum extrapolation of charm contribution to charges
and momentum fraction using Eq. (7), and shown as a function of $M_{\pi}^{2}$
with $a=0$.

Figure 8: Chiral and continuum extrapolation of charm contribution to charges
and momentum fraction using Eq. (7), and shown as a function of lattice
spacing, with $M_{\pi}$ set to $135$ MeV.

Comparison is made with previous results in Tab. 2 for the strange quark and
in Tab. 3 for the charm quark. We find that the axial and tensor charges are
consistent with previous results. Our scalar charge has large uncertainty
because of the smaller statistics on finer lattices. The momentum fraction
(without including mixing correction) is much larger.

strange | $g^{s,R}_{S}$ | $g^{s,R}_{A}$ | $g^{s,R}_{T}$ | $\langle x\rangle_{s}$
---|---|---|---|---
This work (preliminary) | $0.13\sim 0.50$ | $-0.0397(73)$ | $-0.0032(19)$ | $0.109(27)$
PNDME’18 [8, 9] | n/a | $-0.053(8)$ | $-0.0027(16)$ | n/a
$\star$ETMC’19 [10, 11] Nf=4 | $0.454(16)$ | $-0.0458(73)$ | $-0.00268(58)$ | 0.052(12)
$\star$ETMC’19 [10] Nf=2 | $0.371(38)$ | $-0.061(17)$ | $-0.0041(12)$ | n/a
$\chi$QCD’18 [12, 13] | n/a | $-0.035(9)$ | n/a | 0.051(26)
$\star$JLQCD’18 [14] | $0.15(16)$ | $-0.046(28)$ | $-0.012(18)$ | n/a
Table 2: Comparison of our strange quark results with previous calculations of the same quantities. Results labeled with $\star$ are from single lattice spacing. charm | $g^{c,R}_{S}$ | $g^{c,R}_{A}$ | $g^{c,R}_{T}$ | $\langle x\rangle_{c}$
---|---|---|---|---
This work (preliminary) | $0.015(34)$ | $-0.0026(18)$ | $0.0004(13)$ | $0.070(25)$
$\star$ETMC’19 Nf=4 | $0.075(17)$ | $-0.0098(34)$ | $-0.00024(16)$ | $0.019(9)$
$\star$ETMC’19 Nf=2 | $0.059(18)$ | $-0.0065(51)$ | $-0.0060(37)$ | n/a
Table 3: Comparison of our charm quark results with previous calculations of
the same quantities. Results labeled with $\star$ are from single lattice
spacing.

##  4 Conclusion

In this work, we present the calculation of strange and charm contribution to
nucleon charges and momentum fraction. We computed the disconnected diagrams
on 6 lattice ensembles covering 4 lattice spacings and 3 pion masses that
allow us to extrapolate to the physical point. The renormalization constants
are calculated in the RI-sMOM scheme and then matched to $\overline{MS}$
scheme and evolved to $\mu=2$ GeV. We take into account the mixing for scalar
charge between strange and iso-scalar, but have not included the mixing for
momentum fractions. Our results on axial charge and tensor charge are
consistent with previous lattice results. The scalar charge calculation needs
increased statistics on finer lattices. The calculation of momentum fraction
also needs more data and the evaluation of mixing with gluon operator and
other flavors.

## References

  * [1] SciDAC, LHPC, UKQCD Collaboration, R. G. Edwards and B. Joo, The Chroma software system for lattice QCD, Nucl. Phys. Proc. Suppl. 140 (2005) 832, [hep-lat/0409003]. [,832(2004)].
  * [2] R. Babich, J. Brannick, R. Brower, M. Clark, T. Manteuffel, S. McCormick, J. Osborn, and C. Rebbi, Adaptive multigrid algorithm for the lattice Wilson-Dirac operator, Phys. Rev. Lett. 105 (2010) 201602, [arXiv:1005.3043].
  * [3] J. Osborn, R. Babich, J. Brannick, R. Brower, M. Clark, S. Cohen, and C. Rebbi, Multigrid solver for clover fermions, PoS LATTICE2010 (2010) 037, [arXiv:1011.2775].
  * [4] PNDME Collaboration, T. Bhattacharya, V. Cirigliano, S. Cohen, R. Gupta, A. Joseph, H.-W. Lin, and B. Yoon, Iso-vector and Iso-scalar Tensor Charges of the Nucleon from Lattice QCD, Phys. Rev. D 92 (2015), no. 9 094511, [arXiv:1506.06411].
  * [5] HPQCD, UKQCD Collaboration, E. Follana, Q. Mason, C. Davies, K. Hornbostel, G. P. Lepage, J. Shigemitsu, H. Trottier, and K. Wong, Highly improved staggered quarks on the lattice, with applications to charm physics, Phys. Rev. D 75 (2007) 054502, [hep-lat/0610092].
  * [6] MILC Collaboration, A. Bazavov et al., Lattice QCD Ensembles with Four Flavors of Highly Improved Staggered Quarks, Phys. Rev. D 87 (2013), no. 5 054505, [arXiv:1212.4768].
  * [7] A. Hasenfratz and F. Knechtli, Flavor symmetry and the static potential with hypercubic blocking, Phys. Rev. D 64 (2001) 034504, [hep-lat/0103029].
  * [8] R. Gupta, B. Yoon, T. Bhattacharya, V. Cirigliano, Y.-C. Jang, and H.-W. Lin, Flavor diagonal tensor charges of the nucleon from (2+1+1)-flavor lattice QCD, Phys. Rev. D 98 (2018), no. 9 091501, [arXiv:1808.07597].
  * [9] H.-W. Lin, R. Gupta, B. Yoon, Y.-C. Jang, and T. Bhattacharya, Quark contribution to the proton spin from 2+1+1-flavor lattice QCD, Phys. Rev. D 98 (2018), no. 9 094512, [arXiv:1806.10604].
  * [10] C. Alexandrou, S. Bacchio, M. Constantinou, J. Finkenrath, K. Hadjiyiannakou, K. Jansen, G. Koutsou, and A. Vaquero Aviles-Casco, Nucleon axial, tensor, and scalar charges and $\sigma$-terms in lattice QCD, Phys. Rev. D 102 (2020), no. 5 054517, [arXiv:1909.00485].
  * [11] C. Alexandrou, S. Bacchio, M. Constantinou, J. Finkenrath, K. Hadjiyiannakou, K. Jansen, G. Koutsou, H. Panagopoulos, and G. Spanoudes, Complete flavor decomposition of the spin and momentum fraction of the proton using lattice QCD simulations at physical pion mass, Phys. Rev. D 101 (2020), no. 9 094513, [arXiv:2003.08486].
  * [12] Y.-B. Yang, J. Liang, Y.-J. Bi, Y. Chen, T. Draper, K.-F. Liu, and Z. Liu, Proton Mass Decomposition from the QCD Energy Momentum Tensor, Phys. Rev. Lett. 121 (2018), no. 21 212001, [arXiv:1808.08677].
  * [13] J. Liang, Y.-B. Yang, T. Draper, M. Gong, and K.-F. Liu, Quark spins and Anomalous Ward Identity, Phys. Rev. D 98 (2018), no. 7 074505, [arXiv:1806.08366].
  * [14] JLQCD Collaboration, N. Yamanaka, S. Hashimoto, T. Kaneko, and H. Ohki, Nucleon charges with dynamical overlap fermions, Phys. Rev. D 98 (2018), no. 5 054516, [arXiv:1805.10507].