# TQCompressor: improving tensor decomposition methods in neural networks via
permutations

V. Abronin1, A. Naumov1, D. Mazur1, D. Bystrov1, K. Tsarova1,
Ar. Melnikov1, I. Oseledets, S. Dolgov1, R. Brasher1, M. Perelshtein1,†
1 Terra Quantum AG, Kornhausstrasse 25, 9000 St. Gallen, Switzerland
<EMAIL_ADDRESS>
(August 2023)

###### Abstract

We introduce TQCompressor, a novel method for neural network model compression
with improved tensor decompositions. We explore the challenges posed by the
computational and storage demands of pre-trained language models in NLP tasks
and propose a permutation-based enhancement to Kronecker decomposition. This
enhancement makes it possible to reduce loss in model expressivity which is
usually associated with factorization. We demonstrate this method applied to
the GPT-2small [15]. The result of the compression is TQCompressedGPT-2 model,
featuring 81 mln. parameters compared to 124 mln. in the GPT-2small. We make
TQCompressedGPT-2 publicly available. We further enhance the performance of
the TQCompressedGPT-2 through a training strategy involving multi-step
knowledge distillation, using only a 3.1% of the OpenWebText [9].
TQCompressedGPT-2 surpasses DistilGPT-2 [12] and KnGPT-2 [1] in comparative
evaluations, marking an advancement in the efficient and effective deployment
of models in resource-constrained environments.

##  1 Introduction

Advancements in pre-trained language models have significantly impacted the
field of natural language processing (NLP), providing robust solutions for
tasks such as language translation, text summarization, and question
answering. Despite their effectiveness, the substantial computational and
storage demands of these models present a notable challenge. Consequently,
research efforts have increasingly concentrated on model compression
strategies to alleviate these demands. Techniques such as knowledge
distillation [5], post-training quantization [19][18], pruning [19], and
matrix factorization, specifically tensor decomposition [11] [10], are at the
forefront of these endeavors.

Matrix factorization methods offer promising prospects, as demonstrated by the
high compression ratios achieved without significant performance drops in
downstream tasks, such as KroneckerBERT’s [2] success on the GLUE benchmark
[3]. While these methods theoretically reduce FLOPs, the full realization of
these benefits is hindered by current hardware and software limitations, which
are not optimized for sparse structures resulting from matrix factorization.

The future, however, looks bright with potential solutions on the horizon.
Advances in computing architecture, particularly in quantum computing and
adapted GPU designs, are set to bridge this gap. Quantum computing, in
particular, with its proficiency in handling complex computations, aligns well
with the demands of matrix factorization. Research in this area, such as
developments in tensor network computational architectures [21] [22], is
already paving the way. These innovations suggest that the theoretical
advantages of matrix factorization could soon be fully realized,

Some works focus specifically on factorizing weights in Language Models (LMs)
using Kronecker decomposition[14]. They have been applied to such models as
BERT [2] and GPT-2 [1], which are reliable model compression benchmarks.
Kronecker decomposition can drastically reduce the models number of parameters
at the expense of its expressivity. This drawback is then corrected by
training the model on a fraction of the original dataset.

This paper introduces a novel permutation-based enhancement to the Kronecker
decomposition method, reducing the drop in model expressivity commonly
associated with factorization. Our technique, applied to the GPT-2small
model—a standard in compression benchmarks—modifies the embeddings, multi-head
attention (MHA), and feed forward neural network (FFN) modules. We replace
these modules with customized, compressed layers, resulting in a streamlined
model with fewer parameters yet comparable performance.

Our method demonstrates its versatility by being applicable to various neural
network matrix decomposition techniques. By applying permutations to the
weight matrices, we rearrange neuron connections, optimizing the structure for
tensor decomposition without altering the network’s fundamental functionality.
This rearrangement enhances the network’s representational power and accuracy.

We further enhance our model’s performance using a training strategy based on
knowledge distillation [5], which effectively mitigates the effects of
compression while utilizing only a selected small subset of the dataset (3.1%
of OpenWebText). Our approach involves an epoch-by-epoch decomposition
technique during training, which incrementally optimizes the model’s
structure. In comparative evaluations, this method surpasses its direct
competitors and shows superior performance compared to DistillGPT-2 [12].

Our contribution:

  * •

We introduce a new permutation algorithm designed to enhance the efficacy of
matrix factorization methods. This algorithm theoretically reduces the
performance degradation typically associated with these methods, thereby
improving model efficiency without compromising on accuracy.

  * •

Our method’s effectiveness is showcased through its application to the
Kronecker Decomposition method, specifically in the context of compressing the
GPT-2small model. This demonstrates not only a reduction in model size but
also a maintenance of high performance levels, indicative of the algorithm’s
potential applicability across various neural network architectures.

  * •

We have made the trained weights of TQCompressedGPT-2 available along with the
code111The code and model are publicly available:
https://huggingface.co/tq-ag/TQCompressedGPT2
https://github.com/terra-quantum-public/TQCompressedGPT2

Figure 1: This figure illustrates the compression pipeline of a pre-trained
GPT-2small model using the our Decomposition algorithm. The process begins
with the original, uncompressed model on the left. The compression algorithm
is applied in the central part of the diagram, consisting of three main steps:
First, the row-column permutation of weight matrices is performed to improve
the representability of the matrices for decomposition. Second step involves
the Kronecker Decomposition of the permuted weight matrices. Algorithm
performs multiple iterations until the desired level of approximation accuracy
is achieved. The outcome of this process is the compressed model –
TQCompressedGPT-2, shown on the right side of the diagram. This compressed
model retains essential performance characteristics while reducing the overall
number of parameters, thereby making it more efficient for deployment and use
in resource-constrained environments.

##  2 Related works

In this section, we review the related works in the field, with a focus on the
Kronecker decomposition algorithm combined with permutations of weight
matrices, which is the approach employed in our work.

Matrix factorization methods have gained prominence as efficient techniques
for compressing neural networks. These methods aim to approximate the weight
matrices of neural networks with a lower-rank factorization, reducing both the
number of parameters and theoretical computational complexity. Matrix
factorization approaches for neural network compression include, but not
limited to singular value decomposition (SVD) [8], Tensor Train decomposition
[17], Tensor Ring [16] and Kronecker decomposition [14].

Numerous studies have focused on the Kronecker decomposition of weights,
primarily due to their ability to achieve high compression rates while
maintaining, or minimally impacting, model perplexity. The Kronecker
decomposition stands out not only for its compression efficacy but also for
its significant reduction in FLOPs, which is a crucial factor in enhancing
computational efficiency.

A notable milestone in this domain was set by Edalati and Tahaei (2021)[1],
who pioneered the use of Kronecker decomposition for compressing the
GPT-2small model. This model, containing approximately 83 million parameters,
is versatile, having been pre-trained under multiple setups to adapt to
various downstream tasks. Their work demonstrated not just a feasible approach
to compressing a complex model like GPT-2 but also opened avenues for further
research in applying such decompositions more broadly across different model
architectures. This breakthrough is particularly significant given the
increasing demand for efficient, high-performing models in a wide range of
applications, from natural language processing to more computationally
intensive tasks.

While factorized models generally exhibit an increase in perplexity, this can
be effectively countered by fine-tuning the factorized model using a subset of
the original dataset. A prevalent strategy for this fine-tuning process is
Knowledge Distillation, as outlined in [5].

In our research, we also analyze DistilGPT-2 [12], a prominent compressed
variant of the GPT-2small model. DistilGPT-2, with its 82 million parameters,
presents a more compact alternative to the 124 million parameters of the
original GPT-2small, illustrating the potential of model compression
techniques in maintaining performance while reducing size. Notably,
DistilGPT-2 was trained via Knowledge Distillation using the entire
OpenWebText dataset [9], which exemplifies the effectiveness of this approach
in creating efficient yet powerful models. In contrast, our method
demonstrates its efficiency by requiring only small portion of the dataset to
recover expressivity drop, significantly reducing the resources needed while
still achieving comparable performance.

##  3 Methodology

In this section, we expound on our contributions to the field of neural
network compression, detailing both the decomposition method and the
corresponding training procedure. Let $W\in\mathbb{R}^{m\times n}$ denote a
weight matrix within the original neural network architecture.
$P\in\mathbb{R}^{m\times m}$ and $C\in\mathbb{R}^{n\times n}$ represent the
learned permutation matrices for rows and columns, respectively, while
$A\in\mathbb{R}^{m_{1}\times n_{1}}$ and $B\in\mathbb{R}^{m_{2}\times n_{2}}$
describe the Kronecker decomposition matrices for the weight matrix $W$.

###  3.1 TQCompressed decomposition

Building upon the premise that Kronecker products can effectively compress
weights across diverse neural network architectures, we present a novel
decomposition approach. This method focuses on identifying an optimal
permutation of the weight matrix, followed by its Kronecker decomposition
(Figure 1). This problem can be mathematically formulated as:

$\min_{A_{i},B_{i},P,C}|PWC-\sum_{i=1}^{r}A_{i}\otimes B_{i}|_{2}^{2}$ (1)

It is crucial to note that the permutation matrices $P$ and $C$ are binary
matrices, each with a single non-zero element per row, and thus they can be
succinctly expressed as permutation vectors of dimensions $n$ and $m$,
respectively. The introduction of permutations into our method incurs an
additional parameter cost of $n+m$ for each decomposed matrix, resulting in
following total parameter count of decomposed matrix:

$m_{1}n_{1}+m_{2}n_{2}+n+m$ (2)

The underlying intuition of our method is akin to shuffling neurons: permuting
neurons does not alter the network’s structural integrity but renders it more
amenable to factorization, enabling a more efficient approximation of each
layer.

###  3.2 Finding optimal decompositions

Determining optimal decompositions is challenging due to the discrete nature
of permutation matrices. Our iterative algorithm (Algorithm 2) alternates
between optimizing for $P$, $C$ (permutation matrices), and $A$, $B$
(Kronecker decomposition matrices), refining weight approximations.

Optimal permutations are determined in isolation, while the Kronecker matrices
$A$ and $B$ are optimized jointly. Each step concludes by fixing the newly
obtained optima, and the algorithm proceeds iteratively, optimizing one
argument at a time to converge to a solution close to the global optimum.

Initial values for $A$ and $B$ are randomly generated, while $P$ and $C$ begin
as identity matrices. The subsequent sections detail each stage of the
algorithm.

####  3.2.1 Optimal permutation matrices

In this section we state the problem of finding optimal permutation matrices
and explain our solution of it. Explanations are only given for the row-
permutation matrix $P$, but the same reasoning can be applied to the column-
permutation matrix $C$. The problem of finding an optimal permutation matrix
$P$ can be formulated as follows:

$\min_{{P}}\|{P}{W}{C}-{A}\otimes{B}\|^{2}_{2}$ (3)

For simplicity we denote the $WC$ product and the Kronecker product $A\otimes
B$ as $W^{(1)}$ and $W^{(2)}$ respectively:

$\min_{{P}}\|{P}{W^{(1)}}-{W^{(2)}}\|^{2}_{2}$ (4)

In this setting the permutation matrix $P$ can be thought of as a bijective
mapping that pairs rows from $W^{(1)}$ to rows from $W^{(2)}$. Thus the
problem is reduced to finding a one-to-one correspondence between the rows of
$W^{(1)}$ and $W^{(2)}$ that minimizes the total mean squared error between
paired rows. This problem statement is equivalent to the assignment problem
[23] with the following cost matrix Appendix 1. Theorem proof.:

$D_{ij}=\displaystyle\sum_{k}|W^{(1)}_{ik}-W^{(2)}_{jk}|,\
D\in\mathbb{R}^{n\times n}$ (5)

We solve the assignment problem using the Hungarian algorithm [7, 23].

// initialize matrices

$A:=rand\\_matrix(m_{1}\times n_{1})$

$B:=rand\\_matrix(m_{2}\times n_{2})$

$P:=id\\_matrix(m)$

$C:=id\\_matrix(n)$

for k iterations do

$A,B:=kron\\_decomp(PWC)$ //kronecker decomposition using SVD

$Dp:=find\\_D(W\cdot C,A\otimes B)$

$P:=hung\\_alg(Dp)$ //hungarian algorithm

$Dc:=find\\_D(W^{T}\cdot P^{T},(A\otimes B)^{T})^{T}$

$C:=hung\\_alg(Dc)$ //hungarian algorithm

end for

Algorithm 1 Hungarian algorithm for solving 5

####  3.2.2 Optimal Kronecker decomposition

When the optimal values for $P$ and $C$ are found, we can find the optimal
values for $A$ and $B$ by solving the following problem:

$\min_{{A_{i}},{B_{i}}}\|W_{perm}-\sum_{i=1}^{r}A_{i}\otimes B_{i}\|^{2}_{2}$
(6)

Here $W_{perm}$ is the product $WPC$ with fixed optimal values for $P$ and
$C$. This problem can be solved using rank-1 SVD decomposition[8]. We find
that our method outperforms vanilla Kronecker decomposition, while adding an
insignificant overhead to the parameter count and inference costs.

###  3.3 TQCompressed embedding

Embedding layers are fundamental in natural language processing (NLP)
applications, serving as sizable, trainable lookup tables that map discrete
tokens to continuous vector spaces. These layers are mathematically
characterized by an embedding matrix $W_{emb}\in\mathbb{R}^{v\times d}$, where
$v$ is the vocabulary size and $d$ is the embedding dimensionality.

In the architecture of GPT models, two embedding layers are utilized: one for
token representation and another for positional encoding within the sequence.
Our compression technique targets the weight matrices of these embedding
layers along the embedding axis. We apply Kronecker product factorization in
the form of $A^{v\times d/f}\otimes B^{f\times d}$, which allows us to
represent the original embedding matrix $W_{emb}$ as the product of two
smaller matrices $A$ and $B$ with a reduction factor $f$.

###  3.4 TQCompressed transformer

The transformer architecture, pivotal to modern NLP models, is composed of two
main elements: multi-head attention (MHA) and feed-forward network (FFN)
layers. Our compression framework encompasses both these components.

For the MHA layers, the attention mechanism is realized by first projecting
the input through three sets of weights—$W^{Q}$, $W^{K}$, and $W^{V}$—to
create the query, key, and value matrices, respectively. The attention output
$O$ is computed by the equation:

$\displaystyle O$
$\displaystyle=\text{softmax}\left(\frac{QK^{T}}{\sqrt{d_{k}}}\right)V,$ (7)

where $d_{k}$ is the scaling factor, typically the dimensionality of the Key
vectors.

Each head in an MHA layer has its individual set of weight matrices
$W^{Q}_{l}$, $W^{K}_{l}$, and $W^{V}_{l}$, which are conventionally
concatenated across heads:

$\displaystyle W^{{}^{\prime}Q}$
$\displaystyle=\text{Concat}(W_{1}^{Q},\ldots,W_{L}^{Q}),$ (8) $\displaystyle
W^{{}^{\prime}K}$ $\displaystyle=\text{Concat}(W_{1}^{K},\ldots,W_{L}^{K}),$
$\displaystyle W^{{}^{\prime}V}$
$\displaystyle=\text{Concat}(W_{1}^{V},\ldots,W_{L}^{V}),$

where $L$ signifies the total number of attention heads.

In our approach, we perform Kronecker-based decomposition on these
concatenated weights as a whole, rather than individually per head:

$\displaystyle W^{{}^{\prime}Q}$ $\displaystyle\approx P^{Q}(A^{Q}\otimes
B^{Q})C^{Q},$ (9) $\displaystyle W^{{}^{\prime}K}$ $\displaystyle\approx
P^{K}(A^{K}\otimes B^{K})C^{K},$ $\displaystyle W^{{}^{\prime}V}$
$\displaystyle\approx P^{V}(A^{V}\otimes B^{V})C^{V},$

Here, $P^{Q}$, $P^{K}$, and $P^{V}$ are the row permutation matrices, and
$C^{Q}$, $C^{K}$, and $C^{V}$ are the corresponding column permutation
matrices for the query, key, and value weights, respectively. This
factorization allows the model to maintain a high degree of representational
fidelity with a reduced number of parameters.

The output of the MHA layers is then projected using another linear
transformation $W^{O}$, which, along with the subsequent FFN weights $W_{1}$
and $W_{2}$, is also subject to our compression algorithm:

$\displaystyle W^{O}$ $\displaystyle\approx P^{O}(A^{O}\otimes B^{O})C^{O},$
(10) $\displaystyle W_{1}$ $\displaystyle\approx P_{1}(A_{1}\otimes
B_{1})C_{1},$ $\displaystyle W_{2}$ $\displaystyle\approx P_{2}(A_{2}\otimes
B_{2})C_{2},$

Figure. 3 illustrates the GPT-2small weight matrix decomposition scheme,
followed by the knowledge distillation, disclosed in detail in the section
below.

###  3.5 Knowledge distillation

This section delineates the application of knowledge distillation (KD) [5] in
the training process of compressed model. Within the scope of our discussion,
the symbols $S$ (student model) and $T$ (teacher model) are used for
representational clarity.

Knowledge distillation was employed by minimizing the divergence between the
probability distributions of the student’s predictions and the soft targets
provided by the teacher’s predictions. Concretely, we minimized the cross-
entropy loss between the outputs (logits) of the student and teacher models.

The composite loss function amalgamates the standard cross-entropy loss and an
additional term representing the cross-entropy of the logits:

$\mathcal{L}{{}_{\text{total}}}(x,y)=\lambda\mathcal{L}{{}_{\text{Cross-
Entropy}}}(S(x),y)+(1-\lambda)\mathcal{L}_{\text{Logits}}(S(x),T(x))$ (11)

In the loss function $\mathcal{L}_{\text{total}}$, the variable $x$ denotes
the input data that the neural network processes, while $y$ stands for the
corresponding correct labels or targets that the network is intended to
predict. The loss function evaluates the performance of the student model $S$
by comparing its predicted output $S(x)$ against the true labels $y$ and the
teacher model’s predictions $T(x)$, enabling the student to learn and
approximate the teacher’s behavior more effectively. $\lambda$ denotes a
balancing coefficient

Figure 2: Knowledge distillation process for compressing the GPT-2small model
into the TQCompressedGPT-2 variant. A distilled subset of the OpenWebText
corpus serves as the training data for the TQCompressedGPT-2, the student
model, which learns under the guidance of the full-sized GPT-2small, the
teacher model. This process is designed to preserve high performance in the
compressed model by closely matching the teacher’s accuracy.

In a departure from the typical implementation of knowledge distillation, our
approach harnesses only a fractional subset of the dataset—approximately 3.1%
of the original OpenWebText dataset. This strategy significantly enhances
training efficiency.

Figure 2 visualizes the tailored knowledge distillation framework employed in
our study.

###  3.6 Iterative compression with knowledge distillation

The initial application of model compression often results in a non-trivial
performance deficit. To mitigate these effects, previous research has
leveraged knowledge distillation (KD) [2] [1]. Our methodology advances this
technique by sequentially compressing different layers at various stages
throughout the training process. This iterative approach allows the model to
dynamically adjust and recalibrate to the compression of individual layers,
which we have found to be instrumental in achieving more rapid convergence and
improved overall model performance. This iterative process, coupled with
knowledge distillation, enables the rest of the model to adapt seamlessly to
the modifications in the compressed layers.

Algorithm 2 illustrates process for our iterative decomposition with
permutations:

Algorithm 2 Iterative decomposition with permutations and knowledge
distillation

Let $\theta$ denote the model weights

Let $L$ denote the set of all layers in the model

Let $L_{\text{C}}$ be the subset of $L$ containing only layers which we want
to compress

for $i\leftarrow 1$ to num_iter do

Identify the $i$-th subset of layers $L_{i}\subseteq L_{\text{C}}$ for
compression

Find optimal row and column permutations $P$, $C$ for layers in $L_{i}$

Apply Kronecker decomposition to obtain $A_{i}$ and $B_{i}$ for each layer in
$L_{i}$

Update $\theta$ to reflect compressed layers using $P$, $C$, $A_{i}$, and
$B_{i}$

Train the compressed model using knowledge distillation

Update $\theta$ with the learned weights from distillation

end for

return $\theta$

The num_iter parameter specifies the total number of iterations in the
compression and training cycle, with each iteration targeting different model
layers for compression. This granular approach not only refines the
compression process but also aligns the student model more closely with the
teacher’s performance.

Figure 3: Only specific layers undergo a compression process — embedding Layer
(E), feed-forward network (FFN), and multi-head attention layer (MHA). This
process involves applying row-permutation (P) and column-permutation (C)
matrices to the original weight matrices, followed by Kronecker decomposition,
represented by matrices A and B. The classifier outputs from both the original
and the compressed models are then used in a knowledge distillation framework,
where the original GPT-2small model serves as the teacher, and
TQCompressedGPT-2 acts as the student. The distillation process is focused on
aligning the classifier outputs, thereby preserving the performance of the
compressed model relative to its original counterpart.

##  4 Experiments

In our study, we focus on training a compressed version of the GPT-2small
model, which we refer to as TQCompressedGPT-2. This model undergoes training
on approximately 3.1% of the OpenWebText dataset, a significantly smaller
subset compared to the full dataset typically used. We then evaluate our
model’s performance on the Causal Language Modeling task, assessing its
perplexity on benchmark datasets such as Wikitext-2, Wikitext-103[25], and
Lambada[24].

###  4.1 Experimental setup

  * •

Pre-training dataset: We utilize 250K texts from the OpenWebText dataset.

  * •

Downstream task: The primary benchmarking task is causal language modeling
(CLM), with evaluations conducted on the Wikitext-2, Wikitext-103, and Lambada
datasets.

  * •

Decomposition shapes: Our decomposition strategy, detailed in Table 1,
specifically excludes the decomposition of attention layers. This decision was
based on observations that decompressing these layers resulted in a more
significant performance drop compared to other layers. However, by applying a
higher compression ratio to the embedding layer, we managed to achieve a
parameter count comparable to DistilGPT-2 and KnGPT-2.

  * •

Knowledge distillation (KD) scheme: As outlined in Section 3.5, we adopt a
knowledge distillation strategy for training our compressed model.

TQCompressedGPT-2 effectively reduces the parameter count of the original
GPT-2Small from 124 million to 81 million. Our comparative analysis includes
DistillGPT-2 and setups similar to KnGPT-2, which do not utilize permutation
matrices in the decomposition of weight matrices. We maintain the practice of
compressing only odd layers of the GPT-2 architecture, a common approach in
prior works [1]. The specific decomposition shapes vary across different
setups, with further details provided in the results table. All models were
trained using the AdamW optimizer.

One of the key aspects of our approach is the recovery of performance post-
compression. By utilizing just a small fraction of the OpenWebText dataset for
pre-training, we achieve performance levels on par with KnGPT2 and DistilGPT2.
It is important to note that our attempts to replicate the results of the
KnGPT2 paper were unsuccessful. Thus, our baseline for comparison is a GPT-2
model compressed using the Kronecker Decomposition method.

Table 1: Matrix shapes for GPT-$2_{small}$, and compressed version of it - TQCompressedGPT-2 Model | Embedding layer |  Attention layer ($Q,K,V,O$) | FFN layer |  Number of parameters
---|---|---|---|---
GPT-2 small | $50527\times 768$ | $768\times 768$ | $768\times 3072$ | 124 mln.
TQCompressedGPT-2 |  $A:50527\times 192$ $B:1\times 4$ | $768\times 768$ |  $A:768\times 1536$ $B:1\times 2$ | 81 mln.

###  4.2 Results

Our evaluation of language modeling capabilities is centered around measuring
perplexity on the Wikitext-2, Wikitext-103, and Lambada testing subsets. This
assessment also includes models without added permutation matrices to compare
the effectiveness of our approach.

Table 2: Benchmarking results on Wikitext-2, Wikitext-103, and Lambada
datasets (lower scores indicate better performance).
* Denotes the percentage of the OpenWebText dataset used to train the model. Model | Wikitext-103 | Wikitext-2 | Lambada | Dataset* |  Number of parameters
---|---|---|---|---|---
GPT-$2_{small}$ | 29.16 | 24.65 | 45.27 | 100% | 124 mln.
TQCompressedGPT-2 | 40.28 | 32.25 | 64.72 | 3.1% | 81 mln.
KnGPT-2 | 40.97 | 32.81 | 67.62 | 3.1% | 81 mln.
DistilGPT-2 | 44.53 | 36.47 | 75.99 | 100% | 82 mln.

##  5 Conclusion

This paper presented a novel approach to compressing neural networks,
particularly the GPT-2 model, by introducing a permutation-based enhancement
to Kronecker decomposition. Our method stands out for its ability to maintain
the structural integrity and performance of the original model while
significantly reducing its size. This is achieved through a series of
innovative steps: optimal permutations of neuron connections, Kronecker
decomposition, and knowledge distillation.

Our TQCompressedGPT-2 model, a compressed version of GPT-2, demonstrates the
efficacy of our method. We are comparing the performance of TQCompressedGPT-2
to the original GPT-2small, DistillGPT-2 and KnGTP-2. Despite our using only a
tiny fraction of the original dataset for KD, we achieve comparable
performance levels to GPT-2small. Note that both DistillGPT-2 and KnGTP-2 used
far larger portions of the original dataset than we, and yet we have
maintained comparable performance. This not only underscores the potential of
our approach in practical applications but also opens new avenues for further
research. The compressed model, with its fewer parameters, becomes a viable
option for deployment in resource-constrained environments.

Looking forward, the implications of this research are far-reaching. Our
method can potentially be applied to a wide range of neural network
architectures beyond GPT-2, paving the way for more efficient AI models in
various domains. The combination of permutation and tensor decompositions
presents a new paradigm in neural network compression, balancing the trade-off
between model size and performance.

Finally, our research raises intriguing questions about the future of neural
network architecture design. The effectiveness of permutations in improving
the suitability of neural networks for compression suggests that future
architectures could be designed with such optimizations in mind from the
outset. As artificial intelligence continues to evolve, techniques like ours
will be critical in ensuring that advanced models are accessible and practical
for a broader range of applications.

## References

  * [1] Edalati, A., Tahaei, M., Rashid, A., Nia, V. P., Clark, J. J., & Rezagholizadeh, M. (2021). Kronecker decomposition for gpt compression. arXiv preprint arXiv:2110.08152.
  * [2] Tahaei, M. S., Charlaix, E., Nia, V. P., Ghodsi, A., & Rezagholizadeh, M. (2021). Kroneckerbert: Learning Kronecker decomposition for pre-trained language models via Knowledge Distillation. arXiv preprint arXiv:2109.06243.
  * [3] Wang, A., Singh, A., Michael, J., Hill, F., Levy, O., & Bowman, S. R. (2018). GLUE: A multi-task benchmark and analysis platform for natural language understanding. arXiv preprint arXiv:1804.07461.
  * [4] Bacciu, D., & Mandic, D. P. (2020). Tensor decompositions in deep learning. arXiv preprint arXiv:2002.11835.
  * [5] Hinton, G., Vinyals, O., & Dean, J. (2015). Distilling the knowledge in a neural network. arXiv preprint arXiv:1503.02531.
  * [6] Wu, J. N. (2016, February). Compression of fully-connected layer in neural network by Kronecker product. In 2016 Eighth International Conference on Advanced Computational Intelligence (ICACI) (pp. 173-179). IEEE.
  * [7] Bougleux, S., Gaüzere, B., & Brun, L. (2017). A Hungarian algorithm for error-correcting graph matching. In Graph-Based Representations in Pattern Recognition: 11th IAPR-TC-15 International Workshop, GbRPR 2017, Anacapri, Italy, May 16–18, 2017, Proceedings 11 (pp. 118-127). Springer International Publishing.
  * [8] Golub, G. H., & Reinsch, C. (1971). Singular value decomposition and least squares solutions. In Handbook for Automatic Computation: Volume II: Linear Algebra (pp. 134-151). Berlin, Heidelberg: Springer Berlin Heidelberg.
  * [9] Gokaslan, A., & Cohen, V. (2019) OpenWebText corpus. http://Skylion007.github.io/OpenWebTextCorpus.
  * [10] Novikov, A., Podoprikhin, D., Osokin, A., & Vetrov, D. P. (2015). Tensorizing neural networks. Advances in neural information processing systems, 28.
  * [11] Hrinchuk, O., Khrulkov, V., Mirvakhabova, L., Orlova, E., & Oseledets, I. (2019). Tensorized embedding layers for efficient model compression. arXiv preprint arXiv:1901.10787.
  * [12] Sanh, V., Debut, L., Chaumond, J., & Wolf, T. (2019). DistilBERT, a distilled version of BERT: smaller, faster, cheaper and lighter. arXiv preprint arXiv:1910.01108.
  * [13] Thakker, U., Beu, J., Gope, D., Zhou, C., Fedorov, I., Dasika, G., & Mattina, M. (2019). Compressing RNNs for IOT devices by 15-38x using Kronecker products. arXiv preprint arXiv:1906.02876.
  * [14] Atteya, M. J. (2016). Kronecker product with applications. MJ J. Algebra Appl., 1, 1-4.
  * [15] Radford, A., Wu, J., Child, R., Luan, D., Amodei, D., & Sutskever, I. (2019). Language models are unsupervised multitask learners. OpenAI blog, 1(8), 9.
  * [16] Zhao, Q., Zhou, G., Xie, S., Zhang, L., & Cichocki, A. (2016). Tensor Ring decomposition. arXiv preprint arXiv:1606.05535.
  * [17] Oseledets, I. V. (2011). Tensor-Train decomposition. SIAM Journal on Scientific Computing, 33(5), 2295-2317.
  * [18] Xiao, G., Lin, J., Seznec, M., Wu, H., Demouth, J., & Han, S. (2023, July). Smoothquant: Accurate and efficient post-training quantization for large language models. In International Conference on Machine Learning (pp. 38087-38099). PMLR.
  * [19] Dettmers, T., Lewis, M., Belkada, Y., & Zettlemoyer, L. (2022). Llm. int8 (): 8-bit matrix multiplication for transformers at scale. arXiv preprint arXiv:2208.07339.
  * [20] Kwon, W., Kim, S., Mahoney, M. W., Hassoun, J., Keutzer, K., & Gholami, A. (2022). A fast post-training pruning framework for transformers. Advances in Neural Information Processing Systems, 35, 24101-24116.
  * [21] Qu, Z., Deng, L., Wang, B., Chen, H., Lin, J., Liang, L., Li, G., Zhang, Z. & Xie, Y. (2021). Hardware-enabled efficient data processing with Tensor-Train decomposition. IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 41(2), 372-385.
  * [22] Lykov, D., Schutski, R., Galda, A., Vinokur, V., & Alexeev, Y. (2022, September). Tensor network quantum simulator with step-dependent parallelization. In 2022 IEEE International Conference on Quantum Computing and Engineering (QCE) (pp. 582-593). IEEE.
  * [23] Kuhn, H. W. (1955). The Hungarian method for the assignment problem. Naval research logistics quarterly, 2(1-2), 83-97.
  * [24] Paperno, D., Kruszewski, G., Lazaridou, A., Pham, Q. N., Bernardi, R., Pezzelle, S., Baroni, M., Boleda, G. & Fernández, R. (2016). The LAMBADA dataset: Word prediction requiring a broad discourse context. arXiv preprint arXiv:1606.06031.
  * [25] Merity, S., Xiong, C., Bradbury, J., & Socher, R. (2016). Pointer sentinel mixture models. arXiv preprint arXiv:1609.07843.

## Appendix

### Appendix 1. Theorem proof.

We demonstrate the equivalence of minimizing (3) and solving the assignment
problem with matrix (4) as follows:

For any matrix A its Euclidean norm can be presented in such way:

$\|A\|^{2}=tr(A^{T}\cdot A)=tr(A\cdot A^{T})$. With respect to that:

$\displaystyle\|PW^{(1)}-W^{(2)}\|^{2}$
$\displaystyle=tr[(PW^{(1)}-W^{(2)})(W^{(1)T}P^{T}-W^{(2)T})]$
$\displaystyle=tr(PW^{(1)}W^{(1)T}P^{T}-PW^{(1)}W^{(2)T}-W^{(2)}W^{(1)T}P^{T}+W^{(2)}W^{(2)T})$
$\displaystyle=tr(PW^{(1)}W^{(1)T}P^{T})-2tr(PW^{(1)}W^{(2)T})+tr(W^{(2)}W^{(2)T})$
$\displaystyle=\|W^{(1)}\|^{2}-2tr(PW^{(1)}W^{(2)T})+\|W^{(2)}\|$

$\|W^{(1)}\|^{2}+\|W^{(2)}\|^{2}=const$. It is equivalent to minimize the
function that differs by a constant value, so

$(\|PW^{(1)}-W^{(2)}\|^{2}\rightarrow\min)\Leftrightarrow(tr(P\cdot
K)\rightarrow\min)$

where $K:=-2W^{(1)}W^{(2)T}$.

Here, the initial problem (3) is equivalent to minimization of $tr(P\cdot K)$,
where $P$ is a permutation matrix. This is by definition the assignment
problem with the matrix $K$.

This explains the values stored in matrix $K$:

$\displaystyle(W^{(1)}W^{(2)T})_{ij}$ $\displaystyle=w^{(1)}_{i}\cdot
w^{(2)T}_{j}=\sum_{k}W^{(1)}_{ik}\cdot W^{(2)}_{jk}$ $\displaystyle K_{ij}$
$\displaystyle=-2(W^{(1)}W^{(2)T})_{ij}=-2\sum_{k}W^{(1)}_{ik}\cdot
W^{(2)}_{jk}$

The multiplication of $i$-th row of $W^{(1)}$ and $j$-th row of $W^{(2)}$
results in the matrix $K$. Let $D$ be a matrix defined as:

$\displaystyle D_{ij}$ $\displaystyle=\|w^{(1)}_{i}-w^{(2)}_{j}\|^{2}$
$\displaystyle=tr[(w^{(1)}_{i}-w^{(1)}_{j})(w^{(1)T}_{i}-w^{(2)T}_{j})]$
$\displaystyle=tr(w^{(1)}_{i}w^{(1)T}_{i})-2tr(w^{(1)}_{i}w^{(2)T}_{j})+tr(w^{(2)}_{j}w^{(2)T}_{j})$
$\displaystyle=\|w^{(1)}_{i}\|^{2}-2\sum_{k}W^{(1)}_{ik}W^{(2)}_{jk}+\|w^{(2)}_{j}\|^{2}.$

As shown, $K_{ij}=-2\sum_{k}W^{(1)}_{ik}W^{(2)}_{jk}\Rightarrow
D_{ij}=\|w^{(1)}_{i}\|^{2}+K_{ij}+\|w^{(2)}_{j}\|^{2}$.

$tr(P\cdot
D)=\displaystyle\sum_{i=1}^{n}D_{i,p(i)}=\sum_{i=1}^{n}(\|w^{(1)}_{i}\|^{2}+K_{i,p(i)}+\|w^{(2)}_{p(i)}\|^{2})=$
$=\sum_{i=1}^{n}\|w^{(1)}_{i}\|^{2}+tr(P\cdot
K)+\sum_{j=1}^{n}\|w^{(2)}_{j}\|^{2}\displaystyle,$

where $p(i)$ is an index permutation corresponding to permutation matrix $P$.

Again,
$\displaystyle\sum_{i=1}^{n}\|w^{(1)}_{i}\|^{2}+\sum_{j=1}^{n}\|w^{(2)}_{j}\|^{2}=const$,
so minimization of $tr(P\cdot D)$ is equivalent to $tr(P\cdot K)$. Therefore,
minimizing (3) is equivalent to the assignment problem with matrix $D$.


# A discrete relativistic spacetime formalism for 1+1-QED with continuum
limits

Kevissen Sellapillay  Aix-Marseille Université, CPT, Campus de Luminy, case
907, 13288 Marseille, France<EMAIL_ADDRESS>these authors
contributed equally to this work  Pablo Arrighi  Université Paris-Saclay,
CNRS, Laboratoire de recherche en informatique, 91405, Orsay, France and IXXI,
Lyon  these authors contributed equally to this work  Giuseppe Di Molfetta
Aix-Marseille Université, Université de Toulon, CNRS, LIS, 13288, Marseille,
France<EMAIL_ADDRESS>these authors contributed equally to
this work

###### Abstract

We build a quantum cellular automaton (QCA) which coincides with $1+1$ QED on
its known continuum limits. It consists in a circuit of unitary gates driving
the evolution of particles on a one dimensional lattice, and having them
interact with the gauge field on the links. The particles are massive
fermions, and the evolution is exactly $U(1)$ gauge-invariant. We show that,
in the continuous-time discrete-space limit, the QCA converges to the Kogut-
Susskind staggered version of $1+1$ QED. We also show that, in the continuous
spacetime limit and in the free one particle sector, it converges to the Dirac
equation—a strong indication that the model remains accurate in the
relativistic regime.

###### keywords:

Schwinger model, QED, Quantum Cellular Automata, Dirac equation

## Introduction

Quantum physical phenomena can always be modelled classically by means of
matrices and vectors. But, as far as we know, the dimension of these vectors
grows exponentially with the number of particles, making these models
intractable for classical computers. To simulate quantum physical phenomena
efficiently, it seems we have no choice but to harness the laws of quantum
mechanics themselves, as Feynman first suggested [1]. Quantum simulation could
be applied to better understand condensed matter problems [2], simulate
molecules, find ground states of Hamiltonians, or even simulate the dynamics
of quantum field theories (QFT) [3, 4, 5]. It is the latter application that
motivates this paper.

Amongst QFT, gauge theories are of fundamental importance to Physics, as they
capture the fundamental interactions. Some of them have been recast in
discrete space. Lattice QCD [6] is the most famous example as it is
extensively used to obtain theoretical numerical values, to be compared
against experimental values coming out of particle accelerators: this
procedure is partly how physicists are searching for new physics. Simulation
has therefore taken a central role in the scientific method of particle
physics. But these techniques are computationally heavy: finding a way to
simulate lattice gauge theories efficiently and accurately with a quantum
simulation device would be a game changer. Lattice gauge theories are also key
for condensed matter through their application in spin liquids, and for
quantum error correction e.g. via Kitaev’s toric code [7, 8].

The $1+1$ QED, also known as the Schwinger model [9], is a good candidate for
a first step towards the quantum simulation of the dynamics of a gauge theory.
Indeed, it is based on the $U(1)$ gauge group just like $3+1$ QED. It captures
many non trivial physical properties such as a mass gap, fermion confinement
and chiral symmetry breaking [10]. It is exactly solvable in the massless
limit [9]. These features explain why it is often used as a testbed for new
techniques and ideas.

The standard ways to quantum simulate QFT are fundamentally non-relativistic,
as they all begin by expressing the theory in continuous-time discrete-space
Hamiltonian form, using Kogut-Susskind methods [11, 12]. They then map the
matter (fermions) and the gauge field (bosons) onto quantum systems on a
lattice, whose interactions will mimic those of the target Hamiltonian [3].
Sometimes these interactions are implemented as discrete-time products of
quantum gates, but even then these are obtained by approximating the target
Hamiltonian via the Trotter formula, an approximation which remains valid only
in the non-relativistic regime $\Delta_{t}\ll\Delta_{x}$. This approach was
recently realized experimentally on an ion trap architecture [13]. Numerical
techniques exist that come to complement the standard approach, based on
tensors networks. Those use compact, approximate descriptions of quantum
states [14] such as the Density Matrix Renormalization Group (DMRG) [15, 16],
discarding unwanted information about the states as they evolve, so that their
description remain of manageable size—whilst attempting to keep track of the
interesting physical ingredients.

In order to achieve quantum simulation in the relativistic regime
$\Delta_{x}\approx\Delta_{t}$ one must keep space and time on an equal
footing, discretizing both at the same time. This is sometimes referred to as
digital quantum simulation. Digital simulation has indeed been very successful
at describing relativistic particles in different fields [17], but so far it
has not been able to produce simulation scheme for interacting QFT in the
sector of more than 2 particles.

A proposal was made in [18]. From a QCA simulating relativistic Dirac equation
a $U(1)$ gauge-invariant model was obtained. The convergence to the Schwinger
model was not shown however and in fact, using the method presented in this
paper, it can be shown that the Hamiltonian in the discrete space continuous
time limit is not the correct one. We present here the circuit that gives the
correct Hamiltonian in this limit.

In [19, 20] a quantum walk (i.e. the one particle sector of a QCA) was
proposed which unifies non-relativistic analog quantum simulation with
relativistic digital quantum simulation. Just by imposing
$\Delta_{x}=\Delta_{t}^{1-\alpha}$ and tuning the $\alpha$, the operator is
found to have well defined limits lattice fermions both in continuous-time
discrete-space, and the relativistic Dirac equation in the continuous
spacetime limit—a property referred to as plasticity.

The QCA presented in this paper is closely related to these last two models.
It is again based upon a QCA that recovers the relativistic Dirac equation,
extended to become natively discrete gauge-invariant as in [18]. But this
time, the QCA is plastic, allowing us to prove its continuum limit towards
$1+1$ QED, in the regime where $1+1$ QED does have a limit, i.e. the non-
relativistic regime. In other words, we recover the Hamiltonian of the Kogut-
Susskind Schwinger model in the continuous-time discrete-space limit. In the
continuous spacetime limit we show that the QCA yields the Dirac equation in
the free one particle sector, allowing to make the bridge between the non-
relativistic and the relativistic regimes. Altogether, the QCA coincides with
$1+1$ QED on its mathematical continuum limits, whenever these are defined.

The natively discrete digital circuit for staggered Schwinger model we propose
is not seen in the literature. The QCA is staggered which is not usual in the
QCA formalism [21]. One may wonder whether the approach, beyond the quantum
simulation application, could be used to reframe QFT. Indeed, the fact that
the QCA is gauge-invariant by construction, contains explicit relativistic and
non-relativistic limits, is expressible by means of path integrals [22, 23],
suggests that Quantum Computing point of view upon QFT may bring both rigour
and pedagogy to the table—reviving the line of thought initiated by Feynman
with his checkerboard propagator for 1+1 Dirac equation [24].

The paper is organized as follows. We first define the QCA model, that is the
spacetime structure and the gates. Second we show the continuous-time and
discrete-space limit towards the Kogut-Susskind version of the Schwinger
model, by means of the Jordan-Wigner transformation from qubits to fermions.
Third we show that in a continuous spacetime limit, we recover the Dirac
equation for the free one particle sector. Finally, we prove that the model is
gauge-invariant and conclude by giving some perspectives.

## Model

### The Kogut-Susskind staggered version of the Schwinger model

The Schwinger model [9] is a $(1+1)$D model invariant under the $U(1)$ gauge
group. It models spinless electrons and their antiparticles, positrons,
propagating on a 1D lattice and interacting with a $U(1)$ gauge field. We
briefly summarize it by giving its Hamiltonian, which can be written using a
temporal gauge ($A_{0}(x)=0$, and $A(x)=A_{1}(x)$) as :

$H=\int\\!\mathrm{d}x\bigl{(}\psi^{\dagger}(x)\left[\left(i\partial_{x}+igA(x)\right)\sigma_{z}+m\sigma_{x}\right]\psi(x)+\frac{1}{2}E^{2}(x)\bigr{)},$
(1)

where $E(x)$ is the electric field observable at $x$, and $A(x)$ is its
conjugate momentum, meaning

$[A(x),E(y)]=i\delta(x-y).$ (2)

Here $\psi(x)=(\psi_{1}(x),\psi_{2}(x))^{T}$ is a two components fermion field
satisfying

$\\{\psi_{\alpha}(x),\psi^{\dagger}_{\beta}(y)\\}=\delta_{\alpha\beta}\delta(x-y).$
(3)

Figure 1: Kogut-Susskind version of the Schwinger model. The gauge (boson)
field is represented by the links in red. These links are states $\ket{l}$
where $l$ takes value in $\mathbb{Z}$. The operators acting on them are
$L_{p}$ and $e^{\pm i\theta_{p}}$. The matter (fermion) field $\phi_{p}$,
$\phi^{\dagger}_{p}$ are on the nodes in black. Even (odd) sites correspond to
upper (lower) component of a spinor field. $a$ is the lattice spacing.

We now describe the staggered Kogut-Susskind version of the Schwinger model.
This Kogut-Susskind procedure consists in putting fermion fields on the nodes
of an infinite 1D lattice and bosonic gauge fields on the links between them,
as depicted in Fig.1. We can interpret occupied odd sites as electrons and
unoccupied even sites as positrons, therefore particles and antiparticles are
described by a single fermion field [13] (staggered picture). This procedure
allows to give a continuous time, discrete space formulation of a continuous
spacetime model and is also a way to partly resolve the fermion doubling
problem.

The Hamiltonian of this staggered version reads [11, 12]:

$\centering\boxed{H_{S}=\frac{i}{2a}\sum\limits_{p}(\phi^{\dagger}_{p+1}e^{-i\theta_{p}}\phi_{p}-h.c.)+m\sum\limits_{p}(-1)^{p}\phi^{\dagger}_{p}\phi_{p}+\frac{ag^{2}}{2}\sum\limits_{p}L^{2}_{p}},\@add@centering$
(4)

where $a$ is the lattice spacing, $g$ the strength of the interaction (the
charge of the particles), $m$ the mass and $\phi$ the fermion field. To
construct the spinor $\Psi_{\tilde{p}}$ from the fermion fields, we group
fields in pairs $(p,p+1)$, where $p$ is even. In between the fermion fields at
$p$ and $p+1$, there is a gauge field link which takes values in $\mathbb{Z}$.
The operators $e^{\pm i\theta_{p}}$ raises or lowers the value of this link
$[p,p+1]$ such that :

$e^{\pm i\theta_{p}}\ket{l}_{p}=\ket{l\pm 1}_{p}$ (5)

and the electric field is given by $E_{p}=gL_{p}$ where

$L_{p}\ket{l}_{p}=l\ket{l}_{p}.$ (6)

### A quantum cellular automaton for 1+1 QED

The model we propose consists in having one qubit per site $p$ separated by
$\Delta x$ and gauge fields located on the links between each sites, at half
step, modeled by states taking values in $\mathbb{Z}$ [18]. This gauge field
could be experimentally represented by qudits or harmonic oscillators. We
choose the evolution operator of the QCA to be :

$G=\underset{p^{\prime}\text{
even}}{\otimes}W^{*}_{p^{\prime}}e^{-\frac{i}{2}\Delta x\Delta
tg^{2}L^{2}}\underset{p\text{ odd}}{\otimes}W_{p}e^{-\frac{i}{2}\Delta x\Delta
tg^{2}L^{2}}.$ (7)

where

$W_{p}=\begin{pmatrix}I&&0&&0&&0\\\ 0&&e^{-i\zeta}\sin\theta I&&\cos\theta
V_{p+\frac{1}{2}}&&0\\\ 0&&-\cos\theta
V^{\dagger}_{p+\frac{1}{2}}&&e^{i\zeta}\sin\theta I&&0\\\
0&&0&&0&&I\end{pmatrix}.$ (8)

corresponds to the fermions dynamics and the exponential applied before each
$W_{p}$ codes the interaction with the bosons field at the same position. In
particular, as depicted in Fig. 2, the gate $W_{p}$ and the gate $W^{*}_{p}$
are located in the space-time grid in between the qubit at position $(p,p+1)$
respectively for odd and even $p$. Each gate $W$ and its conjugate acts on the
local gauge field. Moreover, in order to avoid the fermion doubling problem we
choose to work with a staggered QCA, in which occupied odd sites are
interpreted as electrons and unoccupied even sites as positrons. This is
consistent with the same interpretation that Kogut and Susskind gave of the
Schwinger model. From a dynamical point of view, the matrix $W_{p}$ can be
interpreted as follows : diagonal terms correspond to staying on the same
site, eventually picking up a phase related to the mass. The non diagonal part
correspond to hopping terms : a right-moving $\ket{1}$ will decrease the gauge
field it passes through, a left-moving one will increase it. The gauge field
operators read:

$V_{p+\frac{1}{2}}\ket{l}_{p+\frac{1}{2}}=\ket{l-1}_{p+\frac{1}{2}},$ (9)

and

$V^{\dagger}_{p+\frac{1}{2}}\ket{l}_{p+\frac{1}{2}}=\ket{l+1}_{p+\frac{1}{2}}.$
(10)

The identities $\mathbb{I}$ means that the gauge link is left invariant. The
operator $L$ in the exponentials of Eq. (7) acts on the states of the gauge
field $\ket{l}_{p+\frac{1}{2}}$ in the following way :

$L\ket{l}_{p+\frac{1}{2}}=l\ket{l}_{p+\frac{1}{2}},$ (11)

for $l\in\mathbb{Z}$.

Figure 2: QCA structure. One application of $G$, the evolution operator,
corresponds to two rows, first a row of $W$ and then of $W^{*}$ gates. Black
wires represent the fermionic fields, green wires represent the gauge field.
The qubits are separated in space by $\Delta x$ and the gates are separated in
time by $\Delta t$.

Now, in order to map the above automaton with the Schwinger model we move to
the second quantization formalism. Let us start with the gate $W_{p}$ (8),
which acts on a pair of qubits $(p,p+1)$ and a gauge link $[p,p+1]$. It is
useful to introduce $\tilde{W}_{p}=I_{<p}\otimes W_{p}\otimes I_{>p}$. Each
$W_{p}$ can be written in terms of single qubit operators
$E_{ij}=\ket{i}\bra{j}$. Using the order
$(0_{p}0_{p+1},0_{p}1_{p+1},1_{p}0_{p+1},1_{p}1_{p+1})$ and omitting
identities on the gauge link for clarity, we have :

$\displaystyle\begin{split}W_{p}=(E_{00})_{p}\otimes(E_{00})_{p+1}+e^{-i\zeta}\sin\theta(E_{00})_{p}\otimes(E_{11})_{p+1}+\\\
e^{i\zeta}\sin\theta(E_{11})_{p}\otimes(E_{00})_{p+1}+\cos\theta
V_{p+\frac{1}{2}}(E_{01})_{p}\otimes(E_{10})_{p+1}\\\ -\cos\theta
V^{\dagger}_{p+\frac{1}{2}}(E_{10})_{p}\otimes(E_{01})_{p+1}+(E_{11})_{p}\otimes(E_{11})_{p+1},\end{split}$
(12)

Thus, we ought to transform qubit operators to fermionic operators, via the
the standard Jordan-Wigner transformations, and finally introduce the
following annihilation and creation operators :

$\displaystyle\begin{split}&\phi_{p}=I_{<p}\otimes(E_{01})_{p}\otimes(\sigma_{z})_{p+1}\otimes(\sigma_{z})_{p+2}...\\\
&\phi^{\dagger}_{p}=I_{<p}\otimes(E_{10})_{p}\otimes(\sigma_{z})_{p+1}\otimes(\sigma_{z})_{p+2}...,\end{split}$
(13)

with $\\{\phi_{p},\phi^{\dagger}_{p^{\prime}}\\}=\delta_{p,p^{\prime}}I$ and
$\\{\phi_{p},\phi_{p^{\prime}}\\}=0$. Notice that the $\sigma_{z}$ to the
right are specifically there to ensure the correct anti-commutation relation
(Supplementary Material 1.), together they form the string operator of the
Jordan-Wigner transformation. Looking at these operators, we see that creating
or annihilating a fermion on a qubit lattice is very much a non-local
operation and we could be worried that the resulting dynamics using such
operators would be non-local and unphysical. However, the terms involving
these operators in the Hamiltonian are quadratic, therefore they become
perfectly local and physical. This can already be seen at the level of the
gates written in this formalism. In fact, for $\tilde{W}_{p}$ we have :

$\displaystyle\begin{split}\tilde{W}_{p}&=\phi_{p}\phi^{\dagger}_{p}\phi_{p+1}\phi^{\dagger}_{p+1}+e^{-i\zeta}\sin\theta\phi_{p}\phi^{\dagger}_{p}\phi^{\dagger}_{p+1}\phi_{p+1}+e^{i\zeta}\sin\theta\phi^{\dagger}_{p}\phi_{p}\phi_{p+1}\phi^{\dagger}_{p+1}\\\
&-\cos\theta V_{p+\frac{1}{2}}\phi_{p}\phi^{\dagger}_{p+1}-\cos\theta
V^{\dagger}_{p+\frac{1}{2}}\phi^{\dagger}_{p}\phi_{p+1}+\phi^{\dagger}_{p}\phi_{p}\phi^{\dagger}_{p+1}\phi_{p+1}.\end{split}$
(14)

Putting together these transformed gates, the global evolution reads:

$\displaystyle\begin{split}G&=\prod\limits_{p^{\prime}\text{
even}}\tilde{W}^{*}_{p^{\prime}}e^{-\frac{i}{2}\Delta x\Delta
tg^{2}L^{2}}\prod\limits_{p\text{ odd}}\tilde{W}_{p}e^{-\frac{i}{2}\Delta
x\Delta tg^{2}L^{2}}.\end{split}$ (15)

## Methods

### Continuous limits

In order to prove that the above staggered QCA, reformulated in terms of
fermionic operators, converges to the Schwinger Hamiltonian, we introduce the
following parametrization :

$\displaystyle\Delta t$ $\displaystyle=$ $\displaystyle\epsilon$
$\displaystyle\Delta x$ $\displaystyle=$ $\displaystyle\epsilon^{1-\alpha}$
$\displaystyle\kappa$ $\displaystyle=$ $\displaystyle\epsilon^{\alpha}$ (16)
$\displaystyle\theta$ $\displaystyle=$ $\displaystyle\arccos(c\kappa)$
$\displaystyle\zeta$ $\displaystyle=$ $\displaystyle
m\frac{(-1)^{\kappa}\epsilon}{\sin(\theta)}$

where the case $\alpha=1$ and $\alpha=0$ correspond respectively to the
continuous time and discrete space limit and the continuous spacetime limit.

In the non relativistic limit, the evolution is continuous in time and
discrete in space. The dynamics is driven by a Hamiltonian, which is recovered
looking at the first order of the global evolution operator of the QCA, as
follows:

$G=e^{-2iH_{QCA}\Delta t}\simeq 1-2i\Delta tH_{QCA}.$ (17)

We then take the limit $\epsilon\rightarrow 0$ on $G$, using the
parametrisation (Continuous limits), for $\alpha=1$. Finally we get :

$\displaystyle\begin{split}\tilde{W}_{p}&\simeq\phi_{p}\phi^{\dagger}_{p}\phi_{p+1}\phi^{\dagger}_{p+1}+(1-im\epsilon)\phi_{p}\phi^{\dagger}_{p}\phi^{\dagger}_{p+1}\phi_{p+1}+(1+im\epsilon)\phi^{\dagger}_{p}\phi_{p}\phi_{p+1}\phi^{\dagger}_{p+1}\\\
&-\epsilon V_{p+\frac{1}{2}}\phi_{p}\phi^{\dagger}_{p+1}-\epsilon
V^{\dagger}_{p+\frac{1}{2}}\phi^{\dagger}_{p}\phi_{p+1}+\phi^{\dagger}_{p}\phi_{p}\phi^{\dagger}_{p+1}\phi_{p+1}\\\
&=1+\epsilon\Big{[}\phi^{\dagger}_{p+1}\phi_{p}V_{p+\frac{1}{2}}-V^{\dagger}_{p+\frac{1}{2}}\phi^{\dagger}_{p}\phi_{p+1}+im(\phi^{\dagger}_{p}\phi_{p}\phi_{p+1}\phi^{\dagger}_{p+1}-\phi_{p}\phi^{\dagger}_{p}\phi^{\dagger}_{p+1}\phi_{p+1})\Big{]}.\end{split}$
(18)

We can rewrite the mass term :

$\displaystyle\begin{split}&\phi^{\dagger}_{p}\phi_{p}\phi_{p+1}\phi^{\dagger}_{p+1}-\phi_{p}\phi^{\dagger}_{p}\phi^{\dagger}_{p+1}\phi_{p+1}=\phi^{\dagger}_{p}\phi_{p}\phi_{p+1}\phi^{\dagger}_{p+1}-(1-\phi^{\dagger}_{p}\phi_{p})\phi^{\dagger}_{p+1}\phi_{p+1}\\\
&\phi^{\dagger}_{p}\phi_{p}(\phi_{p+1}\phi^{\dagger}_{p+1}+\phi^{\dagger}_{p+1}\phi_{p+1})-\phi^{\dagger}_{p+1}\phi_{p+1}=\phi^{\dagger}_{p}\phi_{p}-\phi^{\dagger}_{p+1}\phi_{p+1}.\end{split}$

Finally we have

$\tilde{W}_{p}\simeq
1+\epsilon\Big{[}\phi^{\dagger}_{p+1}\phi_{p}V_{p+\frac{1}{2}}-V^{\dagger}_{p+\frac{1}{2}}\phi^{\dagger}_{p}\phi_{p+1}+im(\phi^{\dagger}_{p}\phi_{p}-\phi^{\dagger}_{p+1}\phi_{p+1})\Big{]},$
(19)

and the interaction with the gauge field can be developped as

$e^{-\frac{i}{2}\Delta x\Delta tg^{2}L^{2}}\simeq 1-\frac{i}{2}\epsilon
g^{2}\sum\limits_{i}L_{i}^{2}.$ (20)

Combining everything, we get for $G$ :

$\displaystyle\begin{split}G&=\prod\limits_{p^{\prime}\text{
even}}\tilde{W}^{*}_{p^{\prime}}e^{-\frac{i}{2}\Delta x\Delta
tg^{2}L^{2}}\prod\limits_{p\text{ odd}}\tilde{W}_{p}e^{-\frac{i}{2}\Delta
x\Delta tg^{2}L^{2}}\\\
&=\prod\limits_{p^{\prime}}(1+\epsilon[\phi^{\dagger}_{2p^{\prime}+1}\phi_{2p^{\prime}}V_{2p^{\prime}+\frac{1}{2}}-\phi^{\dagger}_{2p^{\prime}}\phi_{2p^{\prime}+1}V^{\dagger}_{2p^{\prime}+\frac{1}{2}}-im(\phi^{\dagger}_{2p^{\prime}}\phi_{2p^{\prime}}-\phi^{\dagger}_{2p^{\prime}+1}\phi_{2p^{\prime}+1})])\\\
&(1-\frac{i}{2}\epsilon
g^{2}\sum\limits_{i}L_{i}^{2})\prod\limits_{p}(1+\epsilon[\phi^{\dagger}_{2p+2}\phi_{2p+1}V_{2p+\frac{3}{2}}-\phi^{\dagger}_{2p+1}\phi_{2p+2}V^{\dagger}_{2p+\frac{3}{2}}\\\
&+im(\phi^{\dagger}_{2p+1}\phi_{2p+1}-\phi^{\dagger}_{2p+2}\phi_{2p+2})])(1-\frac{i}{2}\epsilon
g^{2}\sum\limits_{i}L_{i}^{2}).\end{split}$ (21)

A straightforward calculation leads us to the leading order of the series :

$\displaystyle\begin{split}G&\simeq
1+\epsilon\sum\limits_{p}[\phi^{\dagger}_{2p+1}\phi_{2p}V_{2p+\frac{1}{2}}-\phi^{\dagger}_{2p}\phi_{2p+1}V^{\dagger}_{2p+\frac{1}{2}}+\phi^{\dagger}_{2p+2}\phi_{2p+1}V_{2p+\frac{3}{2}}-\phi^{\dagger}_{2p+1}\phi_{2p+2}V^{\dagger}_{2p+\frac{3}{2}}\\\
&-ig^{2}L_{p}^{2}-im(\phi^{\dagger}_{2p}\phi_{2p}-2\phi^{\dagger}_{2p+1}\phi_{2p+1}+\phi^{\dagger}_{2p+2}\phi_{2p+2})].\end{split}$
(22)

Notice that the mass term can be rewritten as :
$-2im(\phi^{\dagger}_{2p}\phi_{2p}-\phi^{\dagger}_{2p+1}\phi_{2p+1})$ Moreover
the four hopping terms simplify because each pair of two terms is separated by
one step :

$\sum\limits_{p}\phi^{\dagger}_{p+1}\phi_{p}V_{p+\frac{1}{2}}-h.c.$ (23)

Finally, we identify $V_{p+\frac{1}{2}}$ to $e^{-i\theta_{p}}$ and
$V^{\dagger}_{p+\frac{1}{2}}$ to $e^{i\theta_{p}}$. Identifying the
Hamiltonian in (22) using $G\simeq 1-2i\epsilon H_{QCA}$, we find the
Hamiltonian of the QCA to be :

$\boxed{H_{QCA}=\sum\limits_{p}[\frac{i}{2}(\phi^{\dagger}_{p+1}\phi_{p}e^{-i\theta_{p}}-h.c)+m(-1)^{p}\phi^{\dagger}_{p}\phi_{p}+\frac{g^{2}}{2}L_{p}^{2}]}.$
(24)

The above one coincides with the Kogut-Susskind Hamiltonian of the Schwinger
model $H_{S}$ with $a=1$. We have thus identified a QCA-based quantum
simulator for a QED toy model, namely a theory for both spinless electrons and
positrons and their interaction with a dynamical gauge field.

We ought to be sure that in the relativistic limit our simulator reproduces
the right dynamics. Here, we give a proof of that in the simplest scenario,
the non-interacting case. Starting from the the one particle sector of the
staggered QCA, we take the $\alpha=0$ relativistic limit, namely for $\Delta
t=\Delta x=\epsilon\rightarrow 0$. (Continuous limits). The local gate which
drives the automaton simplifies as follow

$W^{\prime}_{p}=\begin{pmatrix}1&&0&&0&&0\\\
0&&e^{-i\zeta}\sin\theta&&\cos\theta&&0\\\
0&&-\cos\theta&&e^{i\zeta}\sin\theta&&0\\\ 0&&0&&0&&1\end{pmatrix},$ (25)

and the global evolution operator reads:

$\displaystyle\begin{split}G=\underset{p^{\prime}\text{
even}}{\otimes}W^{*}_{p^{\prime}}\underset{p\text{
odd}}{\otimes}W_{p},\end{split}$ (26)

Let’s start with a general 1 particle state :
$\ket{\Psi(t)}=\sum\limits_{x}\psi(t,x)\ket{1}_{x}$ and separate the one
particle state into pairs of even and odd sites :

$\displaystyle\begin{split}\ket{\Psi(t)}=\sum\limits_{x\text{
even}}\psi^{l}(t,x)\ket{1}_{x}+\psi^{r}(t,x+\epsilon)\ket{1}_{x+\epsilon},\end{split}$
(27)

After one time step evolution of the automaton, the recurrence relations on
the amplitudes $\mathbf{s}\equiv(\psi^{l},\psi^{r})$ which governs the single
fermion, reads:

$\displaystyle\mathbf{s}(t+2\epsilon,x)=$
$\displaystyle\begin{pmatrix}0&-\cos\theta e^{-i\zeta}\sin\theta\\\
0&\cos^{2}\theta\end{pmatrix}\mathbf{s}(t,x-2\epsilon)+\begin{pmatrix}e^{-2i\zeta}\sin^{2}\theta&\cos\theta
e^{i\zeta}\sin\theta\\\ -\cos\theta
e^{-i\zeta}\sin\theta&e^{2i\zeta}\sin^{2}\theta\end{pmatrix}\mathbf{s}(t,x)$
(28) $\displaystyle+\begin{pmatrix}\cos^{2}\theta&0\\\ \cos\theta
e^{i\zeta}\sin\theta&0\end{pmatrix}\mathbf{s}(t,x+2\epsilon).$ (29)

Taking the limit $\epsilon\rightarrow 0$, we find the following differential
equation :

$\partial_{t}\mathbf{s}(t,x)=P\partial_{x}\mathbf{s}(t,x)+Q\mathbf{s}(t,x).$
(30)

where the operators $P$ and $Q$ are represented in the computational basis :

$P=\begin{pmatrix}c^{2}&c\sqrt{1-c^{2}}\\\
c\sqrt{1-c^{2}}&-c^{2}\end{pmatrix}$ (31)

and

$Q=\begin{pmatrix}im\sqrt{1-c^{2}}&-cim\\\
-cim&-im\sqrt{1-c^{2}}\end{pmatrix}.$ (32)

The operator $P$ is self-adoint and its eigenvalues are $\pm c$. Two
eigenvectors associated to these eigenvalues are :

$\begin{split}b_{-}=-\sqrt{\frac{1-c}{2}}b_{0}+\sqrt{\frac{1+c}{2}}b_{1}\hskip
28.45274ptb_{+}=\sqrt{\frac{1+c}{2}}b_{0}+\sqrt{\frac{1-c}{2}}b_{1}\end{split}$
(33)

The family $(b_{-},b_{+})$ forms an orthonormal basis of the two dimensional
spin Hilbert space. Let us now rewrite equation (30) in this new orthonormal
basis. A straightforward computation leads to:

$i\gamma^{0}\partial_{0}\mathbf{\tilde{s}}(t,x)+i\gamma^{1}\partial_{1}\mathbf{\tilde{s}}(t,x)-m\mathbf{\tilde{s}}(t,x)=0,$
(34)

where $\gamma^{0}=\sigma_{x}$, $\partial_{0}=\partial_{t}$,
$\gamma^{1}=\sigma_{x}\sigma_{z}$ and $\partial_{0}=c\partial_{x}$.

### Remark

Let us shortly discuss the last term of operator $W(p)$. In general, when two
creation operators get exchanged, a minus sign is produced so as to respect
the anti-commutation of their creation operators. The abrupt minus appearance
of minus signs makes it harder to compute quantities about fermions—an issue
which is sometimes referred to as the ‘sign problem’. Here, when $\Delta
t\approx\Delta x$, we must indeed put a $-1$ in the bottom right coefficient
of the gate, since two creation operators are crossing during the lapse of one
$W$ gate, as was shown in detail [18]. However, when $\Delta t\ll\Delta x$ we
should not. Physically, this is because the two fermions now hardly have time
to cross in the lapse of one gate. Mathematically, this shows through the fact
that placing a $-1$ at this position forbids the development of the gate
around identity, ruining any effort to obtain a non-relativistic continuous-
time discrete-space limit towards the Kogut-Susskind Hamiltonian, or any other
Hamiltonian for that matter. To get the best of both worlds, we use scaling
factor $e^{i(\frac{\Delta t}{\Delta x})^{2}\pi}=e^{i\epsilon^{2\alpha}\pi}$,
making the coefficient go to $1$ in the non-relativistic parametrization
($\alpha=1$) and to $-1$ in the relativistic parametrization ($\alpha=0$).
This is compatible with unitarity, plasticity, and fermionic computation.

$W^{\prime\prime}_{p}=\begin{pmatrix}1&&0&&0&&0\\\ 0&&\sin\theta&&-\cos\theta
V_{p+\frac{1}{2}}&&0\\\ 0&&\cos\theta
V^{\dagger}_{p+\frac{1}{2}}&&\sin\theta&&0\\\
0&&0&&0&&e^{i\epsilon^{2\alpha}\pi}\end{pmatrix}.$ (35)

It is the choice that yields the Kogut-Susskind Hamiltonian in the non-
relativistic regime. In the relativistic regime, however, $\alpha=0$ and so at
order 0, we see that :

$\displaystyle\begin{split}\tilde{W}_{p}&\simeq\phi_{p}\phi^{\dagger}_{p}\phi_{p+1}\phi^{\dagger}_{p+1}+\phi_{p}\phi^{\dagger}_{p}\phi^{\dagger}_{p+1}\phi_{p+1}+\phi^{\dagger}_{p}\phi_{p}\phi_{p+1}\phi^{\dagger}_{p+1}-\phi^{\dagger}_{p}\phi_{p}\phi^{\dagger}_{p+1}\phi_{p+1}\\\
&=1-2\phi^{\dagger}_{p}\phi_{p}\phi^{\dagger}_{p+1}\phi_{p+1}.\end{split}$
(36)

Since this is a control-Z on the $2$ qubits, the order $0$ of $G$ cannot be
the identity. As expected, we cannot recover a many-body interacting
Hamiltonian in this regime.

## The staggered QCA is gauge invariant

In this section, we show that our QCA-based quantum simulator is gauge
invariant at finite scale. We define the gauge transformation in Fig.3, where
the gates are defined as follows :

$P_{\varphi}=\underset{p}{\otimes}P_{\varphi(p)}=\underset{p}{\otimes}\Big{(}T_{\varphi(p)}\otimes
R_{\varphi(p)}\otimes T_{-\varphi(p)}\Big{)},$ (37)

where $R_{\varphi(p)}$ acts on the qubit at site $p$ such that :

$\displaystyle R_{\varphi(p)}:$ $\displaystyle\ket{0}\rightarrow\ket{0}$ (38)
$\displaystyle\ket{1}\rightarrow e^{i\varphi(p)}\ket{1},$ (39)

and $T_{\varphi(p)}$ acts on the gauge field states such that :

$T_{\varphi(p)}\ket{l}=e^{il\varphi(p)}\ket{l}.$ (40)

From definition (37), we see that a gauge field at a given site
$p+\frac{1}{2}$ will be acted upon twice, once by $-\varphi(p)$ and once by
$+\varphi(p+1)$. Is our QCA invariant under the above gauge transformation, or
namely:

$P_{\varphi}G\stackrel{{\scriptstyle?}}{{=}}GP_{\varphi}.$ (41)

Figure 3: The gauge transformation. Black dots are the fermions sites. Green
dots in between sites represent the gauge field values. At a single site $p$,
the gauge transformation will apply a gate $R_{\varphi(p)}$ that gives a phase
when the qubit there is in state $\ket{1}$ and apply gates $T_{\varphi(p)}$ on
the left and right gauge field points that produce phases according to the
values of the gauge field. This transformation is applied on every sites,
therefore a single gauge field is acted upon twice.

Equation (41) means that whatever the field $\varphi(p)$ is, the evolution
operator should give the same result if we apply it after the gauge
transformation or before the gauge transformation, in other words, the
dynamics should not be changed by this gauge field, the latter should have no
physical consequences. However if we impose this $U(1)$ local phase,
invariance is only possible if we introduce some interaction between qubits
and the gauge fields , which is done through the $V$ and $V^{\dagger}$
operators, that will compensate the missing pieces of phase appearing during
the particle movement, as done in [18]. To show the gauge invariance (41), we
start from a ket $\ket{1}$ at a given site (odd or even), with only gauge
field values given in the rest of the space. We compute explicitly the results
after applying first $P_{\varphi}$ then $G$ and the other way around, and show
that the result are exactly the same, therefore the gauge transformation
$P_{\varphi}$ has no physical consequences on the dynamics $G$. The model is
thus gauge invariant (Supplementary Material 2.).

## Conclusion

We have described a quantum cellular automaton (QCA) that simulates $1+1$ QED.
It consists in a lattice of qubits encoding whether a fermion is present at a
given site. These interact with the gauge field that lives on the links
between those sites. The QCA was shown to coincide with the Kogut-Susskind
version of the Schwinger model in the continuous-time discrete-space limit,
and with the Dirac equation in the continuous spacetime limit in the one-
particle sector. We go from one limit to the other just by imposing
$\Delta_{x}=\Delta_{t}^{1-\alpha}$ and tuning the $\alpha$.

We still cannot ascertain the QCA recovers the Schwinger model in the
interactive regime of the continuous spacetime limit, but then again it is not
even clear that the lattice QFT has such a limit in the first place. The QCA
coincides with the Schwinger model wherever it has a known, mathematically
defined limit.

It coincides in the story it tells; of fermions propagating relativistically
updating the gauge field, which in turn simply triggers a phase—thereby
turning on the interaction.

Finally, it coincides in terms of its construction: the QCA retains the
fundamental $U(1)$ gauge-invariance of QED, even in discrete spacetime. This
gauge-invariance construction was originally proposed in [18] for QCA, and in
[25] for Reversible CA. Both drew inspirations from gauge-invariant Quantum
Walks [26, 17, 27]. The question of extending these constructions to non-
abelian gauge theories has been treated for Quantum Walks [28, 29] and for
Reversible CA [30]. Lifting this to obtain $U(N)$ gauge-invariant QCA is no
doubt one of the next steps lying ahead towards digital quantum simulation
schemes for Yang Mills theories, such as QCD. An observation was made in [31]
that the path integral of such a theory can be written as a real-time transfer
operator which is itself a finite-depth local quantum circuit.

Another obvious next step is the extension to $2+1$ and $3+1$ dimensions. Of
particular concern is the fact that we relied upon the Jordan Wigner mapping
to encode fermions into qubits. This transformation is known to suffer
locality issues [32, 33, 34], which did not affect us because our terms were
quadratic, cancelling out all non-local effects. In further dimensions we may
not be so lucky however, and will have to rely on alternative transformations
[35, 36, 37]. A possible route of investigation could be to use formulations
of QED in 2+1D or 3+1 dimensions where the Hilbert space dimension is reduced
by either approximating the gauge group [38] or by going in a rotating frame
[39] which decouples the matter to the gauge field and keeps only local
constraints on the latter.

The problem of preparing the ground state of such QCA is puzzling as was
pointed out in [18]. We do not solve the problem but by recovering a
Hamiltonian in the continuous-time discrete-space limit, we make the problem
well-defined.

QCA are closely related to path integrals [22]. Lately a formalism was
developed for dealing with interactions in a perturbative manner [23], within
the QCA framework—which would be an interesting application here. Another vast
topic for exploration is trying to understand the link between renormalization
theory, and the way the parameters of the unitary gates must be made to vary
with the lattice as we take our limits.

The plastic Quantum Walk [19] upon which this plastic QCA is built is
deformable to the point that a curved spacetime limit can been obtained just
by putting a spacetime dependence in the parameter $c$. A fair question to ask
is whether allowing for the same spacetime dependence here, would yield the
Schwinger model on a curved background [40].

## Acknowledgements

We benefited from discussions with Cédric Bény, Nathanaël Éon, Terry Farrelly,
Alberto Verga. This publication was made possible through the support of the
ID# 61466 grant from the John Templeton Foundation, as part of the “The
Quantum Information Structure of Spacetime (QISS)” Project (qiss.fr). The
opinions expressed in this publication are those of the author(s) and do not
necessarily reflect the views of the John Templeton Foundation.

## References

  * [1] Feynman, R. P.  Simulating physics with computers.  _International Journal of Theoretical Physics_ 21, 467–488, DOI: 10.1007/BF02650179 (1982).
  * [2] Lloyd, S.  Universal Quantum Simulators.  _Science_ 273, 1073–1078, DOI: 10.1126/science.273.5278.1073 (1996).  Publisher: American Association for the Advancement of Science _eprint: https://science.sciencemag.org/content/273/5278/1073.full.pdf.
  * [3] Bañuls, M. C. _et al._ Simulating lattice gauge theories within quantum technologies.  _The European Physical Journal D_ 74, 165, DOI: 10.1140/epjd/e2020-100571-8 (2020).
  * [4] Preskill, J.  Simulating quantum field theory with a quantum computer.  In _Proceedings of The 36th Annual International Symposium on Lattice Field Theory — PoS (LATTICE 2018)_ , 024, DOI: 10.22323/1.334.0024 (Sissa Medialab, Michigan State University, East Lansing, Michigan, USA., 2019).
  * [5] Jordan, S. P., Lee, K. S. & Preskill, J.  Quantum algorithms for quantum field theories.  _Science_ 336, 1130–1133 (2012).
  * [6] Knechtli, F., Günther, M. & Peardon, M.  _Lattice Quantum Chromodynamics: Practical Essentials_.  SpringerBriefs in Physics (Springer Netherlands, 2017).
  * [7] Savary, L. & Balents, L.  Quantum spin liquids: a review.  _Reports on Progress in Physics_ 80, 016502, DOI: 10.1088/0034-4885/80/1/016502 (2016).  Publisher: IOP Publishing.
  * [8] Kitaev, A. Y.  Fault-tolerant quantum computation by anyons.  _Annals of Physics_ 303, 2–30, DOI: 10.1016/S0003-4916(02)00018-0 (2003).
  * [9] Schwinger, J.  Gauge Invariance and Mass. II.  _Physical Review_ 128, 2425–2429, DOI: 10.1103/PhysRev.128.2425 (1962).
  * [10] Melnikov, K. & Weinstein, M.  Lattice Schwinger model: Confinement, anomalies, chiral fermions, and all that.  _PHYSICAL REVIEW D_ 14 (2000).
  * [11] Kogut, J. & Susskind, L.  Hamiltonian formulation of Wilson’s lattice gauge theories.  _Physical Review D_ 11, 395–408, DOI: 10.1103/PhysRevD.11.395 (1975).
  * [12] Banks, T., Susskind, L. & Kogut, J.  Strong-coupling calculations of lattice gauge theories: (1 + 1)-dimensional exercises.  _Physical Review D_ 13, 1043–1053, DOI: 10.1103/PhysRevD.13.1043 (1976).
  * [13] Martinez, E. A. _et al._ Real-time dynamics of lattice gauge theories with a few-qubit quantum computer.  _Nature_ 534, 516–519, DOI: 10.1038/nature18318 (2016).  Number: 7608 Publisher: Nature Publishing Group.
  * [14] Orús, R.  A practical introduction to tensor networks: Matrix product states and projected entangled pair states.  _Annals of Physics_ 349, 117–158, DOI: 10.1016/j.aop.2014.06.013 (2014).
  * [15] Byrnes, T. M. R., Sriganesh, P., Bursill, R. J. & Hamer, C. J.  Density matrix renormalization group approach to the massive Schwinger model.  _Physical Review D_ 66, 013002, DOI: 10.1103/PhysRevD.66.013002 (2002).
  * [16] Zapp, K. & Orús, R.  Tensor network simulation of QED on infinite lattices: Learning from $(1+1)\text{ }\mathrm{d}$, and prospects for $(2+1)\text{ }\mathrm{d}$.  _Physical Review D_ 95, 114508, DOI: 10.1103/PhysRevD.95.114508 (2017).  Publisher: American Physical Society.
  * [17] Di Molfetta, G., Brachet, M. & Debbasch, F.  Quantum walks in artificial electric and gravitational fields.  _Physica A: Statistical Mechanics and its Applications_ 397, 157–168 (2014).
  * [18] Arrighi, P., Bény, C. & Farrelly, T.  A quantum cellular automaton for one-dimensional QED.  _Quantum Information Processing_ 19, 88, DOI: 10.1007/s11128-019-2555-4 (2020).
  * [19] Di Molfetta, G. & Arrighi, P.  A quantum walk with both a continuous-time limit and a continuous-spacetime limit.  _Quantum Information Processing_ 19, 47, DOI: 10.1007/s11128-019-2549-2 (2019).
  * [20] Manighalam, M. & Di Molfetta, G.  Continuous time limit of the dtqw in 2d+ 1 and plasticity.  _Quantum Information Processing_ 20, 1–24 (2021).
  * [21] Arrighi, P.  An overview of quantum cellular automata.  _Natural Computing_ 18, 885–899, DOI: 10.1007/s11047-019-09762-6 (2019).
  * [22] D’Ariano, G. M., Mosco, N., Perinotti, P. & Tosini, A.  Path-sum solution of the Weyl Quantum Walk in 3+1 dimensions.  _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 375, 20160394, DOI: 10.1098/rsta.2016.0394 (2017).  ArXiv: 1705.08552.
  * [23] Bisio, A., Mosco, N. & Perinotti, P.  Scattering and perturbation theory for discrete-time dynamics.  _arXiv:1912.09768 [quant-ph]_ (2019).  ArXiv: 1912.09768.
  * [24] Feynman, R. P., Hibbs, A. R. & Styer, D. F.  _Quantum Mechanics and Path Integrals_ (Courier Corporation, 2010).
  * [25] Arrighi, P., Di Molfetta, G. & Eon, N.  A gauge-invariant reversible cellular automaton.  In _International Workshop on Cellular Automata and Discrete Complex Systems_ , 1–12 (Springer, 2018).
  * [26] Di Molfetta, G. & Debbasch, F.  Discrete-time quantum walks: Continuous limit and symmetries.  _Journal of Mathematical Physics_ 53, 123302 (2012).
  * [27] Márquez-Martín, I., Arnault, P., Di Molfetta, G. & Pérez, A.  Electromagnetic lattice gauge invariance in two-dimensional discrete-time quantum walks.  _Physical Review A_ 98, 032333 (2018).
  * [28] Arnault, P., Di Molfetta, G., Brachet, M. & Debbasch, F.  Quantum walks and non-Abelian discrete gauge theory.  _Physical Review A_ 94, 012335, DOI: 10.1103/PhysRevA.94.012335 (2016).
  * [29] Di Molfetta, G. & Pérez, A.  Quantum walks as simulators of neutrino oscillations in a vacuum and matter.  _New Journal of Physics_ 18, 103038 (2016).
  * [30] Arrighi, P., Molfetta, G. D. & Eon, N.  Non-abelian gauge-invariant cellular automata.  In _Theory and Practice of Natural Computing \- 8th International Conference, TPNC 2019, Kingston, ON, Canada, December 9-11, 2019, Proceedings_ , 211–221, DOI: 10.1007/978-3-030-34500-6_15 (2019).
  * [31] Farrelly, T. & Streich, J.  Discretizing quantum field theories for quantum simulation.  _arXiv:2002.02643 [hep-lat, physics:hep-th, physics:quant-ph]_ (2020).  ArXiv: 2002.02643.
  * [32] Huerta, L. & Zanelli, J.  Bose-Fermi transformation in three-dimensional space.  _Physical Review Letters_ 71, 3622–3624, DOI: 10.1103/PhysRevLett.71.3622 (1993).
  * [33] D’Ariano, G. M., Manessi, F., Perinotti, P. & Tosini, A.  Fermionic computation is non-local tomographic and violates monogamy of entanglement.  _EPL (Europhysics Letters)_ 107, 20009, DOI: 10.1209/0295-5075/107/20009 (2014).  ArXiv: 1307.7902.
  * [34] D’Ariano, G. M., Manessi, F., Perinotti, P. & Tosini, A.  The Feynman problem and Fermionic entanglement: Fermionic theory versus qubit theory.  _International Journal of Modern Physics A_ 29, 1430025, DOI: 10.1142/S0217751X14300257 (2014).  ArXiv: 1403.2674.
  * [35] Farrelly, T.  Insights from Quantum Information into Fundamental Physics.  _arXiv:1708.08897 [quant-ph]_ (2017).  ArXiv: 1708.08897.
  * [36] Brun, T. A. & Mlodinow, L.  Quantum cellular automata and quantum field theory in two spatial dimensions.  _Physical Review A_ 102, 062222, DOI: 10.1103/PhysRevA.102.062222 (2020).  ArXiv: 2010.09104.
  * [37] Mlodinow, L. & Brun, T. A.  Fermionic and bosonic quantum field theories from quantum cellular automata in three spatial dimensions.  _arXiv:2011.05597 [quant-ph]_ (2020).  ArXiv: 2011.05597.
  * [38] Haase, J. F. _et al._ A resource efficient approach for quantum and classical simulations of gauge theories in particle physics.  _Quantum_ 5, 393, DOI: 10.22331/q-2021-02-04-393 (2021).  ArXiv: 2006.14160.
  * [39] Bender, J. & Zohar, E.  Gauge redundancy-free formulation of compact QED with dynamical matter for quantum and classical computations.  _Physical Review D_ 102, 114517, DOI: 10.1103/PhysRevD.102.114517 (2020).
  * [40] Parker, L. & Toms, D.  _Quantum Field Theory in Curved Spacetime: Quantized Fields and Gravity_ (Cambridge University Press, 2009).