# Transforming Collections of Pauli Operators into Equivalent Collections of
Pauli Operators over Minimal Registers

Lane G. Gunderman<EMAIL_ADDRESS>The Institute for Quantum
Computing, University of Waterloo, Waterloo, Ontario, N2L 3G1, Canada
Department of Physics and Astronomy, University of Waterloo, Waterloo,
Ontario, N2L 3G1, Canada

(August 27, 2024)

###### Abstract

Transformations which convert between Fermionic modes and qubit operations
have become a ubiquitous tool in quantum algorithms for simulating systems.
Similarly, collections of Pauli operators might be obtained from solutions of
non-local games and satisfiability problems. Drawing on ideas from
entanglement-assisted quantum error-correcting codes and quantum convolutional
codes, we prove the obtainable lower-bound for the number of qubits needed to
represent such Pauli operations which are equivalent and provide a procedure
for determining such a set of minimal register Pauli operations.

††preprint: APS/123-QED

The initial idea for a quantum computer was argued as a tool for simulating
quantum systems [1]. This problem was advanced by Lloyd proposing a method for
actually simulating quantum systems by approximating it with products of
operators [2]. Unfortunately such products have very poor error-scaling unless
higher order product series are used such as those of Suzuki [3, 4].
Additionally, this seminal work by Lloyd did not specify what sort of
operations would be used in the product formula. The simplest building block
operations for quantum devices are the Pauli operators, so if the system can
be simulated using these operations the feasibility of utilizing some sort of
near-term advantage would be improved, providing a useful algorithm for NISQ
(noisy intermediate-scale quantum) devices [5]. A variety of transformations
have been devised for transforming Hamiltonians for quantum systems, such as
molecules or field-theories, into Pauli operators with associated weights [6,
7, 8, 9]. Well controlled qubits with long coherence time are a precious
resource, especially in more near-term quantum devices. These all perform
reasonably well in this regard, however, this work solves one open question
related to these transformations: for a particular collection of Pauli
operators obtained from these transformations, or from some other setting,
what is the minimal number of qubits needed in order to perform such
operations? Our answer is found from quantum error-correction techniques.
Additionally, our solution finds the obtainable lower-bound and also provides
an efficient methodology for transforming the original collection of Pauli
operators into an equivalent collection using the minimal number of qubits.
Note that this work has more general applications–wherever Pauli operators
appear (from quantum complexity theory to non-local games to fault-tolerance
to communication protocols)–but we focus here on more near-term algorithms
related to simulating quantum systems.

##  I Background Definitions

Formally, quantum computers operate by performing unitary operators. The
rudimentary operations for qubits are the Pauli operators, which act as
follows on the computational basis states $\\{|0\rangle,|1\rangle\\}$:

$X|j\rangle=|j\oplus 1\rangle,\quad Z|j\rangle=(-1)^{j}|j\rangle,\quad
j\in\\{0,1\\}.$ (1)

There is one final non-trivial operator, $Y$, which we take as the product of
$X$ and $Z$ 111While $XZ=-ZX$ and $Y=iXZ$, for the sake of our arguments it
suffices to note that $Y$ is just some scalar times the product. The $n$-fold
tensor product of Pauli operators is indicated by $\mathbb{P}^{n}$.

Although the Pauli operators form a basis, working with these non-commuting
operators is rather cumbersome and so we utilize the symplectic representation
for the Pauli operators to transform the tool kit from those of non-commuting
operator theory to those of linear algebra over $\mathbb{Z}_{2}$ [11, 12].

######  Definition 1.

The binary symplectic representation of a Pauli operator is the homomorphic
map, $\phi$, which maps $\mathbb{P}^{n}\mapsto\mathbb{Z}_{2}^{2n}$. It takes
the $i$-th register in the Pauli and sets the $i$-th position in the binary
vector to the power of the $X$ operator at that position and the $i+n$-th
position to the power of the $Z$ operator. In short, this performs:

$\phi\left(\bigotimes_{t=1}^{n}X^{a_{t}}Z^{b_{t}}\right)=\left(\bigoplus_{t=1}^{n}a_{t}\right)\bigoplus\left(\bigoplus_{t=1}^{n}b_{t}\right).$
(2)

This mapping observes:

$\phi(s_{i}\circ s_{j})=\phi(s_{i})\oplus\phi(s_{j})$ (3)

where $s_{i}$ and $s_{j}$ are two Paulis in $\mathbb{P}^{n}$ and $\oplus$ is
componentwise vector addition mod two.

######  Definition 2.

The commutator between two Pauli operators
$\phi(p_{1})=(\vec{a}_{1}|\vec{b}_{1})$ and
$\phi(p_{2})=(\vec{a}_{2}|\vec{b}_{2})$ is computed as:

$\phi(p_{1})\odot\phi(p_{2})=\vec{a}_{1}\cdot\vec{b}_{2}-\vec{b}_{2}\cdot\vec{a}_{2}\mod
2$ (4)

This is $0$ if $p_{1}$ and $p_{2}$ commute, and $1$ if not (they anti-
commute).

######  Definition 3.

For a collection of Pauli operators $\mathcal{P}$ we denote by
$rank_{\phi}(\mathcal{P})$ the number of compositionally independent
generators for $\mathcal{P}$, we call this the $\phi$-rank of $\mathcal{P}$.

The $\phi$-rank of $\mathcal{P}$ is simply the rank of the symplectic
representation of $\mathcal{P}$. Since there will be a second rank of
importance in this work, we have made this separate distinction.

##  II Auxiliary Results

In this section we begin by providing an example illustrating the result of
this work, then state the main result, following it by proving smaller results
in this section before showing the main result in the next.

######  Example 4.

Consider the following pair of operators $\\{XX,IZ\\}$ as a motivating
example. These two operators anti-commute and the rank of their symplectic
form is $2$, however, a shorter set of Pauli operators which also obey these
relations is $X,Z$. So the number of physical registers needed can be reduced
for these operators from $2$ to just $1$. This work considers this problem and
shows the optimal reduction in the number of registers for collections of
Pauli operators so that the same commutation relations are preserved.

######  Definition 5.

Supplementary Pauli operators are a collection of Pauli operators,
$\mathcal{P}^{\cup}$, associated with each Pauli in $\mathcal{P}$ such that
$\mathcal{P}^{\cup}\otimes_{r}\mathcal{P}$ has all operators commuting, where
$\otimes_{r}$ is the tensor product of element $i$ from $\mathcal{P}^{\cup}$
with element $i$ from $\mathcal{P}$.

######  Example 6.

Consider again the prior example, then we have $\mathcal{P}=\\{XX,IZ\\}$ and
$\mathcal{P}^{\cup}=\\{X,Z\\}$ with
$\mathcal{P}^{\cup}\otimes_{r}\mathcal{P}=\\{XXX,ZIZ\\}$.

######  Remark 7.

Note that
$\phi(\mathcal{P}_{i})\odot\phi(\mathcal{P}_{j})=\phi(\mathcal{P}^{\cup}_{i})\odot\phi(\mathcal{P}^{\cup}_{j})$
since we are working in the qubit, binary, case. Generally
$\phi(\mathcal{P}_{i})\odot\phi(\mathcal{P}_{j})+\phi(\mathcal{P}^{\cup}_{i})\odot\phi(\mathcal{P}^{\cup}_{j})=0$,
with $+$ being standard addition mod the local-dimension value [12].

Given this, we will focus on finding a set of $\mathcal{P}^{\cup}$ that has a
minimal number of registers. This problem is exactly that which entanglement-
assisted quantum error-correcting codes needed to solve–finding a set of
supplemental Paulis to act on pre-shared entanglement to resolve non-commuting
generators for the codes [13, 14, 15]. This same problem also appears in
determining the number of registers that must be carried over between frames
in quantum convolutional codes [16]. Fortunately, this problem has been
studied and been solved, although such solutions do not provide a method for
determining these supplementary Paulis [17, 16]. For these works the primary
concern was generating Pauli operators that with their supplementary Pauli
operators pairwise commuted, whereas here we focus on ensuring these
collections have the same commutation relationships. The universality of this
result from EAQECC and quantum convolutional codes is brought forth in this
work, showing that the utility of this method extends beyond error-correction
but can also play a central role in improving some simulation based quantum
algorithms and any place where stabilizers are a powerful tool. Additionally,
as it provides a tight lower-bound for the Paulis obtained from a fermionic to
qubit transformation, it shows that further reduction is not possible, so
other avenues must be considered, such as clique formation [18, 19, 20, 21,
22, 23] or randomized Trotterization [24, 25] or some other method which maps
the Hamiltonian to Pauli operators.

Our primary result is then the following Theorem, which we will prove
following auxiliary Lemmas and associated Definitions for these Lemmas.

######  Theorem 8.

Given a collection of Pauli operators with associated weights on $n$
registers, equivalent Pauli operators with the same $\phi$-rank can be
efficiently determined which obtain the minimal number of registers needed.

In order to reduce this problem, we may consider just a basis subset from the
collection of Pauli operators $\mathcal{P}$:

######  Lemma 9.

It suffices to find generators for the smallest subgroup that contains all of
the given Pauli operators, $\mathcal{P}$.

###### Proof.

Let $\mathcal{P}\subseteq S$ and $s$ be a set of generators for $S$ such that
$s\subseteq\mathcal{P}$. Since the homomorphism $\phi$ carries composition as
addition and the symplectic product is linear in the $\phi$ representations
for Pauli operators, the elements of $\mathcal{P}$ that are not in $s$ can be
generated from composition combinations of elements in $s$. Thus to fully
specify $\mathcal{P}$, $s$ suffices–the weights associated with each operator
must be tracked separately, however, these determine the full set of
commutation relations that the collection of Pauli operators must satisfy. ∎

######  Remark 10.

This result holds for a collection of supplementary Pauli operators,
$\mathcal{P}^{\cup}$.

Note that since we wish to also satisfy the commutation relations also for
compositions of the generators, we require
$rank_{\phi}(\mathcal{P})=rank_{\phi}(\mathcal{P}^{\cup})$.

######  Definition 11.

The commutation matrix for a collection of Pauli operators $\mathcal{P}$,
$M(\mathcal{P})$, has entries
$[M(\mathcal{P})]_{ij}=\phi(\mathcal{P}_{i})\odot\phi(\mathcal{P}_{j})$ for
basis members of $\mathcal{P}$.

From Remark 7, we know that $M(\mathcal{P})=M(\mathcal{P}^{\cup})$.

######  Lemma 12.

The commutator matrix’s rank for a collection of Paulis $\mathcal{P}$,
$rank(M(\mathcal{P}))$, is invariant under the choice of basis elements.

###### Proof.

Let us take some note on how composing generators alters the commutator matrix
$M$. The composition of generator $i$ onto generator $j$, written as $i\circ
j$, results in the simultaneous addition of row $i$ in $M$ to row $j$, as well
as the addition of column $i$ to column $j$ in $M$. This is evident from the
homomorphism rule of composition being equivalent to addition of vectors as
well as addition in the commutators.

Generally, the simultaneous addition row and column $i$ to $j$ is equal to the
row addition, then the column addition, then the addition of $[M]_{ii}$ to
$[M]_{jj}$. These commutativity matrices are anti-symmetric (including the
qudit case) and have diagonals that are all zero, as each operator trivially
commutes with itself, which ensures that the recipient diagonal entry is left
as zero. Then for these matrices the simultaneous addition of rows and columns
is equal to row addition followed by column addition. These are elementary
matrix operations so they will preserve the rank. ∎

##  III Determining a valid set of minimal register operators

As the commutation relations between a collection of Pauli operators
determines how the members of this weighted collection interact with each
other, finding any collection of Pauli operators which satisfy the same
commutation relations will perform the same operations as the original
collection. This means that by processing the collection of Paulis before
running the experiment we may reduce the number of registers needed in the
experiment, then upon completing said experiment covert the results back in
terms of the original operators and in turn the original Hamiltonian.
Alternatively, it states that one can use whichever method one desires to
determine a collection of Pauli operators which satisfy the commutativity
requirements then apply this method to compress those operators into ones
which utilize the fewest qubits.

In a prior work it was found that the size of the minimal-memory encoder for a
quantum convolutional code $S$ was given by $dim(M(S))-\frac{1}{2}rank(M(S))$
222Note that the rank is computed over mod $2$, as all operations are
performed over this field. [16], which is equivalent to the number of
registers needed in our supplementary Paulis. This formula can be derived by
breaking the commutation matrix into its isotropic portion (similar to a null
space) and symplectic portions (a direct sum of rank $2$ matrices). While this
formula provided the number of qubits needed, it did not provide any method
for obtaining the memory operators themselves. Algorithm 1 from [16] provides
a procedure for shortening the size of the memory operators, effectively
removing any additional memory operators needed that may have been included
due to composition of delayed generators, however that is not relevant in this
work. While a useful tool, these results do not provide a methodology for
generating a set of valid memory operators. The following Lemma provides such
a constructive method, which is of particular use here. We will use this Lemma
in the proof of Theorem 8.

######  Lemma 13.

There is a deterministic, efficient procedure for determining a valid set of
supplementary Pauli operators $\mathcal{P}^{\cup}$ using the minimal number of
registers for a given collection of Pauli operators $\mathcal{P}$ with
commutation matrix $M$.

###### Proof.

Let $dim(M)=d$. From Lemma 12 the composition of generators is equivalent to
conjugation by elementary row additions. Upon composing generators we also
retain the anti-symmetric nature of $M$. Begin with the first nonzero column
in $M$ and add rows such that this column is given by $0100\ldots$. The first
row will also be in this form. Repeat this procedure now with the second
column, which will have a nonzero first entry. Continue in this way,
performing Gaussian elimination, until the matrix is, up to permutations,
given by:

$\tilde{D}:=\left(\bigoplus_{i=0}^{dim(M)-rank(M)}[0]\right)\bigoplus\left(\bigoplus_{i=0}^{\frac{1}{2}rank(M)}\begin{bmatrix}0&1\\\
1&0\end{bmatrix}\right).$ (5)

Let $L^{-1}$ be the set of elementary operations performed to transform $M$ to
$\tilde{D}$–so that $M=L\tilde{D}L^{T}$. A set of Pauli operators with
commutativity matrix $\tilde{D}$ is given by
$\\{Z_{i}\\}_{i=1}^{dim(M)-rank(M)}$ and
$\\{X_{i},Z_{i}\\}_{i=dim(M)-rank(M)+1}^{dim(M)-\frac{1}{2}rank(M)}$. To
recover a set of Paulis satisfying $M$ with the minimal number of registers,
$\mathcal{P}^{\cup}$, we may simply compose these aforementioned generators’
$\phi$ representations with $L$. ∎

Putting together the auxiliary results from the prior section and Lemma 13, we
have proven Theorem 8.

######  Corollary 14.

This result holds for qudits as well.

Note that this Lemma in essence re-proves the result from [17, 13] in a
modified language, circumventing the need to appeal to symplectic vector
spaces and the Symplectic Graham-Schmidt Orthogonalization Procedure (SGSOP),
instead working with nearly standard Gaussian elimination. This also means
that this is likely the computationally fastest routine that can exist for
determining the $\mathcal{P}^{\cup}$ as well as the equivalent problems of
finding the priorly shared entanglement in EAQECC and the memory operators in
quantum convolutional codes. Additionally this proof provides a methodology
for determining a satisfying set of Pauli operators, a piece that was lacking
in prior works, as the reduction occurs through repeated row additions
(compositions of the generators to be found) but does not require any Clifford
operations.

To close this section we provide an example illustrating this result:

######  Example 15.

We will take:

$\mathcal{P}=\\{ZYZZXZXYXI,IIXYYZIYYY,\\\ IYXIXIXYZY,YZZIXZZXYI,ZZZYIXYXXZ,\\\
XZZIXIIXIZ,XXIYYIIYIX,IXIXIYIIYI\\}.$ (6)

There are $8$ operators on $10$ qubits. From the $\phi$ representation, we
know that we will need between $4$ and $8$ registers (inclusive) since those
are the minimal rank while using both the $X$ and $Z$ components of the
representation and the maximal rank when only needing one component (all
operators commute). For these Paulis the commutator matrix is given by:

$M(\mathcal{P})=\begin{bmatrix}0&0&0&1&1&1&0&0\\\ 0&0&0&1&1&0&1&0\\\
0&0&0&1&0&0&1&0\\\ 1&1&1&0&0&1&0&0\\\ 1&1&0&0&0&1&0&0\\\ 1&0&0&1&1&0&0&1\\\
0&1&1&0&0&0&0&1\\\ 0&0&0&0&0&1&1&0\end{bmatrix}$ (7)

This matrix has (binary) rank $6$, which from our formula says that we should
need exactly $5$ registers to represent this set of commutation relations.

We perform the following operations to put the matrix into our pseudodiagonal
form: add $5$ to $4$, then add $4$ to $7$, then add $7$ to $6$, then add $5$
to $6$, then add $8$ to $2$, then add $6$ to $1$, then add $6$ to $2$:

$\begin{bmatrix}0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0\\\ 0&0&0&1&0&0&0&0\\\
0&0&1&0&0&0&0&0\\\ 0&0&0&0&0&1&0&0\\\ 0&0&0&0&1&0&0&0\\\ 0&0&0&0&0&0&0&1\\\
0&0&0&0&0&0&1&0\end{bmatrix}.$ (8)

This matrix is now in the reduced form and can be satisfied by
$\\{Z_{1},Z_{2},X_{3},Z_{3},X_{4},Z_{4},X_{5},Z_{5}\\}$. Then matrix $L$ which
we conjugated by to obtain this form is given by the inverse of the product of
these elementary row additions, providing:

$L=\begin{bmatrix}1&0&0&0&0&1&0&0\\\ 0&1&0&0&0&1&0&1\\\ 0&0&1&0&0&0&0&0\\\
0&0&0&1&1&0&0&0\\\ 0&0&0&0&1&0&0&0\\\ 0&0&0&0&1&1&1&1\\\ 0&0&0&1&0&0&1&0\\\
0&0&0&0&0&0&0&1\end{bmatrix}.$ (9)

This satisfies $M(\mathcal{P})=L\tilde{D}L^{T}$, so then a minimal register
collection of Pauli operators can be given by:

$\mathcal{P}^{\cup}=\\{Z_{1}Z_{4},Z_{2}Z_{4}Z_{5},X_{3},Z_{3}X_{4},X_{4},Y_{4}Y_{5},Z_{3}X_{5},Z_{5}\\}$
(10)

These Pauli operators on $5$ registers, with weight at most $3$, satisfy the
same commutation matrix $M$ for $\mathcal{P}$, so these form a minimal
register set of equivalent Pauli operators.

##  IV Conclusion

Pauli operators form the most fundamental building blocks in quantum
computations, error-correction, and algorithms. In this work we have couched
the discussion of the possible uses of this result in terms of near-term
algorithms, however, it can be of rather general use, since knowing the
minimal number of registers to represent a collection of commutation relations
could be useful for a variety of settings. This work poses the question of
what the minimal number of registers needed to represent a particular
commutation relation is, then proceeds to solve this problem. This work proves
an optimal register compression algorithm for a given set of commutation
relations. In particular, this work imports a central idea from entanglement-
assisted quantum error-correction and quantum convolutional codes into the
broader framework of Pauli operators for any purpose. From the subset of
examples considered thus far on small molecules, it seems that the Bravyi-
Kitaev transformation achieves this optimal register usage 333Andrew Jena,
private communications[6]. Is this always true and if so, can this be proven?
Another future direction to carry this work would be to see whether this
method, with appropriate extra considerations, can be used to prove lower-
bounds on the weights of Pauli operators needed in order to obtain a
collection of commutation relations, and how close a deterministic procedure
can get to this minimal weight. Another important avenue to consider would be
how the results are impacted by connectivity constraints of a physical device.

## Acknowledgments

We thank Andrew Jena for helpful discussions, David Cory for feedback, and
Mark Wilde for convincing the author to look more closely at EAQECC and
quantum convolutional codes.

## References

  * [1] Richard P Feynman.  Simulating physics with computers.  In Feynman and computation, pages 133–153. CRC Press, 2018.
  * [2] Seth Lloyd.  Universal quantum simulators.  Science, 273(5278):1073–1078, 1996.
  * [3] Masuo Suzuki.  Fractal decomposition of exponential operators with applications to many-body theories and monte carlo simulations.  Physics Letters A, 146(6):319–323, 1990.
  * [4] Masuo Suzuki.  General theory of fractal path integrals with applications to many-body theories and statistical physics.  Journal of Mathematical Physics, 32(2):400–407, 1991.
  * [5] John Preskill.  Quantum computing in the nisq era and beyond.  Quantum, 2:79, 2018.
  * [6] Sergey B Bravyi and Alexei Yu Kitaev.  Fermionic quantum computation.  Annals of Physics, 298(1):210–226, 2002.
  * [7] Andrew Tranter, Peter J Love, Florian Mintert, and Peter V Coveney.  A comparison of the bravyi–kitaev and jordan–wigner transformations for the quantum simulation of quantum chemistry.  Journal of chemical theory and computation, 14(11):5617–5630, 2018\.
  * [8] Jarrod R McClean, Ryan Babbush, Peter J Love, and Alán Aspuru-Guzik.  Exploiting locality in quantum computation for quantum chemistry.  The journal of physical chemistry letters, 5(24):4368–4380, 2014\.
  * [9] Jarrod R McClean, Nicholas C Rubin, Kevin J Sung, Ian D Kivlichan, Xavier Bonet-Monroig, Yudong Cao, Chengyu Dai, E Schuyler Fried, Craig Gidney, Brendan Gimby, et al.  Openfermion: the electronic structure package for quantum computers.  Quantum Science and Technology, 5(3):034014, 2020.
  * [10] While $XZ=-ZX$ and $Y=iXZ$, for the sake of our arguments it suffices to note that $Y$ is just some scalar times the product.
  * [11] Michael A Nielsen and Isaac Chuang.  Quantum computation and quantum information, 2002.
  * [12] Avanti Ketkar, Andreas Klappenecker, Santosh Kumar, and Pradeep Kiran Sarvepalli.  Nonbinary stabilizer codes over finite fields.  IEEE Transactions on Information Theory, 52(11):4892–4914, 2006\.
  * [13] Todd Brun, Igor Devetak, and Min-Hsiu Hsieh.  Correcting quantum errors with entanglement.  science, 314(5798):436–439, 2006.
  * [14] Min-Hsiu Hsieh.  Entanglement-assisted coding theory.  arXiv preprint arXiv:0807.2080, 2008.
  * [15] Charles H Bennett, David P DiVincenzo, John A Smolin, and William K Wootters.  Mixed-state entanglement and quantum error correction.  Physical Review A, 54(5):3824, 1996.
  * [16] Monireh Houshmand, Saied Hosseini-Khayat, and Mark M Wilde.  Minimal-memory, noncatastrophic, polynomial-depth quantum convolutional encoders.  IEEE transactions on information theory, 59(2):1198–1210, 2012\.
  * [17] Mark M Wilde and Todd A Brun.  Optimal entanglement formulas for entanglement-assisted quantum coding.  Physical Review A, 77(6):064302, 2008.
  * [18] Vladyslav Verteletskyi, Tzu-Ching Yen, and Artur F Izmaylov.  Measurement optimization in the variational quantum eigensolver using a minimum clique cover.  The Journal of chemical physics, 152(12):124114, 2020.
  * [19] Andrew Jena, Scott Genin, and Michele Mosca.  Pauli partitioning with respect to gate sets.  arXiv preprint arXiv:1907.07859, 2019.
  * [20] Ariel Shlosberg, Andrew J Jena, Priyanka Mukhopadhyay, Jan F Haase, Felix Leditzky, and Luca Dellantonio.  Adaptive estimation of quantum observables.  arXiv preprint arXiv:2110.15339, 2021.
  * [21] Danny Paulson, Luca Dellantonio, Jan F Haase, Alessio Celi, Angus Kan, Andrew Jena, Christian Kokail, Rick van Bijnen, Karl Jansen, Peter Zoller, et al.  Simulating 2d effects in lattice gauge theories on a quantum computer.  PRX Quantum, 2(3):030334, 2021.
  * [22] Ewout van den Berg and Kristan Temme.  Circuit optimization of hamiltonian simulation by simultaneous diagonalization of pauli clusters.  Quantum, 4:322, 2020.
  * [23] Pranav Gokhale, Olivia Angiuli, Yongshan Ding, Kaiwen Gui, Teague Tomesh, Martin Suchara, Margaret Martonosi, and Frederic T Chong.  Minimizing state preparations in variational quantum eigensolver by partitioning into commuting families.  arXiv preprint arXiv:1907.13623, 2019.
  * [24] Earl Campbell.  Random compiler for fast hamiltonian simulation.  Physical review letters, 123(7):070503, 2019.
  * [25] Andrew M Childs, Aaron Ostrander, and Yuan Su.  Faster quantum simulation by randomization.  Quantum, 3:182, 2019.
  * [26] Note that the rank is computed over mod $2$, as all operations are performed over this field.
  * [27] Andrew Jena, private communications.


# A tensor formalism for multilayer network centrality measures using the
Einstein product

Smahane El-Halouy11footnotemark: 1<EMAIL_ADDRESS>Silvia
Noschese22footnotemark: 2<EMAIL_ADDRESS>Lothar
Reichel33footnotemark: 3<EMAIL_ADDRESS>Department of Mathematical
Sciences, Kent State University, Kent, OH 44242, USA, and Laboratory LAMAI,
Faculty of Sciences and Technologies, Cadi Ayyad University, Marrakech,
Morocco.  Dipartimento di Matematica “Guido Castelnuovo”, SAPIENZA Università
di Roma, P.le A. Moro, 2, I-00185 Roma, Italy.  Department of Mathematical
Sciences, Kent State University, Kent, OH 44242, USA.

###### Abstract

Complex systems that consist of different kinds of entities that interact in
different ways can be modeled by multilayer networks. This paper uses the
tensor formalism with the Einstein tensor product to model this type of
networks. Several centrality measures, that are well known for single-layer
networks, are extended to multilayer networks using tensors and their
properties are investigated. In particular, subgraph centrality based on the
exponential and resolvent of a tensor are considered. Krylov subspace methods
are introduced for computing approximations of different measures for large
multilayer networks.

###### keywords:

multilayer networks , centrality measures , adjacency tensor , tensor
functions , Einstein product , Krylov subspace method

###### MSC:

[2010] 05C50 , 15A18 , 65F15

##  1 Introduction

A network is a set of objects that are connected to each other in some
fashion. Mathematically, a _single-layer network_ is represented by a graph
$G=\\{V,E\\}$, where the elements of the set $V=\\{v_{i}\\}_{i=1}^{n}$,
referred to as vertices or nodes, represent the objects, and the elements of
the set $E=V\times V$, designated as edges, represent the connections between
the nodes. We denote an edge from node $v_{i}$ to node $v_{j}$ by
$v_{i}\rightarrow v_{j}$.

Some real world examples require the modeling of more than one kind of nodes
or of more than one type of edges. This holds, for instance, for the
transportation network in a country when considering different means of
transportation. The train and bus routes are different types of connections
and should in some models be represented by different kinds of edges.
Moreover, train and bus stations may make up nodes with diverse properties.
The connections between a train station and an adjacent bus station give rise
to yet another kind of edges connecting different kinds of nodes, along which
travelers typically walk. This kind of objects and connections can be modeled
by _multilayer networks_ , which emphasize different kinds or connections,
known as layers, between possibly different kinds of elements of a network.
Each layer is represented by a single graph that contains the elements, or
some of the elements, of the network and the connections between them in this
level. Edges connecting nodes from different layers model the interactions
between different levels. Therefore, the nodes in a multilayer network require
two indices, e.g., $v_{i}^{\ell}$, where the superscript $\ell$ denotes the
layer, and the subscript $i$ determines the node in this layer. The set
$V_{L}=V\times L$, where $V$ is the set of nodes in each layer and $L$ is the
set of layers, represents all possible combinations of node-layers, and the
set of edges $E=V_{L}\times V_{L}$ corresponds to pairs of possible
combinations of nodes and layers. The special case when the set of nodes is
the same in all layers, and edges that connect nodes in different layers are
only allowed between a node and its copy in another layer, is known as a
_multiplex network_. A nice recent paper by Bergermann and Stoll [7] studies
multiplex networks and generalized matrix function-based centrality measures
to thid kind of networks. The authors use supra-adjacency matrices to
represent multiplex networks. Recently, a global measure of communicability in
a multiple network, computed by means of the Perron root, and the right and
left Perron vectors of the supra-adjacency matrix associated with the network
was introduced in [15].

The multilayer networks mentioned above can be generalized to allow for the
modeling of more than one property of the nodes in each layer. Therefore, we
define a sequence $\\{L_{j}\\}_{j=1}^{d}$ of sets of layers with $d$ being the
number of different properties that one would like to model; $L_{j}$ is the
set of layers for the property $j$. Then $V_{L}=V\times
L_{1}\times\ldots\times L_{d}$ and the nodes are identified as
$v_{i}^{\ell_{1},\ldots,\ell_{d}}$. For more details on this kind of
generalization, we refer to [11, 29] and the references therein, where general
frameworks for multilayer network are discussed together with their
mathematical formulation.

Single-layer networks are often represented by an adjacency matrix, which is
helpful for extracting information about the network by evaluating matrix
functions or by computing certain eigenvectors. For instance, Estrada and
Higham [19] describe how the matrix exponential and resolvent can be used to
determine how easy it is to communicate between nodes in a single-layer
network, and which nodes are the most important ones; see also Estrada [17]
and references therein. For multilayer networks, we use tensors, i.e., a
multidimensional generalization of matrices, and represent the network by a
$2(d+1)$-order adjacency tensor $\mathcal{A}$ of size
$(|V|\times|L_{1}|\times\ldots\times|L_{d}|)\times(|V|\times|L_{1}|\times\ldots\times|L_{d}|)$,
where $|V|$ denotes the number of nodes in each layer, and $|L_{j}|$, with
$j=1,2,\ldots,d$, designates the number of layers. The entry
$\mathcal{A}(i,\ell_{1},\ldots,\ell_{d},j,k_{1},\ldots,k_{d})$ of the tensor
$\mathcal{A}$ for an _unweighted_ network with layers $L_{j}$,
$j=1,2,\ldots,d$ is one if there is an edge
$v_{i}^{\ell_{1},\ldots,\ell_{d}}\rightarrow v_{j}^{k_{1},\ldots,k_{d}}$;
otherwise the tensor entry is zero. For a _weighted_ network the entry entry
is replaced by a real, generally positive, number. In a directed network some
of the edges represent “one-way streets”.

Adjacency tensors allow us to capture the structure and complexity of
relationships between the nodes of a multilayer network. We are interested in
investigating and generalizing some centrality measures that are well
established for single-layer networks to multilayer networks by using the
tensor formalism and applying tensor tools, such as the tensor product and
tensor functions. Several centrality measures have been studied for multilayer
networks using the tensor formalism in [12]. Eigenvector multicentrality has
been investigated for multilayer networks via a tensor-based framework in
[33]; see also [10]. In addition to generalizing centrality measures that are
commonly used for single-layer networks to multilayer networks, we describe
practical and efficient ways to compute these measures by using Krylov
subspace methods.

This paper is organized as follows. Tensor notation, definitions, and
properties used throughout this paper are described in Section 2. Section 3
discusses the extension of matrix functions to tensor functions by using the
tensor Einstein product. We define centrality measures for multilayer networks
based on the tensor representation and tensor functions. Section 4 describes
Krylov subspace methods based on the tensor format using the Einstein product,
and discusses their application to the approximation of tensor functions.
Section 5 presents a few computed examples and Section 6 contains concluding
remarks.

##  2 Preliminaries

This section presents notation and properties of tensors that will be used
throughout this paper. We start with a generalization of the matrix-matrix
product to tensors that is referred to as the Einstein product.

######  Definition 1 (Einstein Product)

Let $\mathcal{A}\in\mathbb{R}^{I_{1}\times\ldots\times I_{N}\times
J_{1}\times\ldots\times J_{M}}$ and
$\mathcal{B}\in\mathbb{R}^{J_{1}\times\ldots\times J_{M}\times
K_{1}\times\ldots\times K_{L}}$ be tensors of orders $N+M$ and $M+L$,
respectively. The product
$\mathcal{C}=\mathcal{A}\ast_{M}\mathcal{B}\in\mathbb{R}^{I_{1}\times\ldots\times
I_{N}\times K_{1}\times\ldots\times K_{L}}$ of the tensors $\mathcal{A}$ and
$\mathcal{B}$ is a tensor of order $N+L$ with entries

$\mathcal{C}_{i_{1},\ldots,i_{N},k_{1},\ldots,k_{L}}=\sum_{j_{1},\ldots,j_{M}}\mathcal{A}_{i_{1},\ldots,i_{N},j_{1},\ldots,j_{M}}\mathcal{B}_{j_{1},\ldots,j_{M},k_{1},\ldots,k_{L}}.$

It is commonly referred to as the Einstein product; see [8, 9, 14]. The
subscript $M$ in $\ast_{M}$ shows the last and first $M$ dimensions of
$\mathcal{A}$ and $\mathcal{B}$, respectively, over which the sum is
evaluated.

The identity tensor
$\mathcal{I}=[\mathcal{I}_{i_{1},\ldots,i_{N},j_{1},\ldots,j_{N}}]\in\mathbb{R}^{I_{1}\times\ldots\times
I_{N}\times I_{1}\times\ldots\times I_{N}}$ under the Einstein product has the
entries

$\mathcal{I}_{i_{1},\ldots,i_{N},j_{1},\ldots,j_{N}}=\left\\{\begin{array}[]{cl}1,&~{}{\rm
if}~{}i_{k}=j_{k},~{}{\rm for}~{}k=1,2,\ldots,N,\\\ 0,&~{}{\rm
otherwise}.\end{array}\right.$

###### Remark 1

A tensor of even order $\mathcal{A}\in\mathbb{R}^{I_{1}\times\cdots\times
I_{N}\times J_{1}\times\cdots\times J_{N}}$ is said to be _square_ if the
first set of dimensions equals the second set, i.e., if $I_{k}=J_{k}$ for
$k=1,2,\ldots,N$; see, e.g., [31]. The adjacency tensor of a multilayer
network is of even order and square, however, for some computations tensors of
different orders are required. This is possible when using the Einstein
product by choosing the number of dimensions over which we carry out the
summation.

###### Remark 2

The transpose of a tensor $\mathcal{A}\in\mathbb{R}^{I_{1}\times\ldots\times
I_{N}\times J_{1}\times\ldots\times J_{M}}$ is a tensor
$\mathcal{B}\in\mathbb{R}^{J_{1}\times\ldots\times J_{M}\times
I_{1}\times\ldots\times I_{N}}$ such that
$\mathcal{A}_{i_{1},\ldots,i_{N},j_{1},\ldots,j_{M}}=\mathcal{B}_{j_{1},\ldots,j_{M},i_{1},\ldots,i_{N}}$;
see, e.g., [31]. A well-known tensor-matrix product is the so-called $n$-mode
product [28], which is the same as the Einstein product when summing over the
last mode of the tensor. In other words, we have

$\mathcal{A}\times_{N+1}A^{T}=\mathcal{A}\ast_{1}A,$

where $\mathcal{A}\in\mathbb{R}^{I_{1}\times\ldots\times I_{N}\times J}$,
$A\in\mathbb{R}^{J\times I}$, and $\times_{N+1}$ denotes the $n$-mode product
over mode $N+1$.

We can reorganize the entries of a tensor in different ways to obtain a $2D$
array, i.e., a matrix. This transformation is known as _matricization_ or
_flattening_. We adopt the following definition:

######  Definition 2 (Tensor flattening)

Let $\mathcal{A}\in\mathbb{R}^{I_{1}\times\ldots\times I_{N}\times
J_{1}\times\ldots\times J_{M}}$ be a tensor of order $N+M$. The elements of
the matrix $A\in\mathbb{R}^{I_{1}\ldots I_{N}\times J_{1}\ldots J_{M}}$
obtained by flattening the tensor $\mathcal{A}$ are given by

$A_{i,j}=\mathcal{A}_{i_{1},\ldots,i_{N},j_{1},\ldots,j_{M}},$

where

$\displaystyle i$ $\displaystyle=$ $\displaystyle
i_{1}+\sum_{p=2}^{N}(i_{p}-1)\prod_{q=1}^{p-1}I_{q},$ $\displaystyle j$
$\displaystyle=$ $\displaystyle
j_{1}+\sum_{p=2}^{M}(j_{p}-1)\prod_{q=1}^{p-1}J_{q}.$

Here and below the indices $i_{p}$ and $j_{p}$ live in their domains, i.e.,
$1\leq i_{p}\leq I_{p}$ for $1\leq p\leq N$, and $1\leq j_{q}\leq J_{p}$ for
$1\leq q\leq M$. We define $A={\rm mat}(\mathcal{A})$ and $\mathcal{A}={\rm
mat}^{-1}(A)$.

###### Remark 3

For a multiplex network, ${\rm mat}(\mathcal{A})$ is the supra-adjacency
matrix defined in [7].

###### Proposition 1

Let $\mathcal{A}\in\mathbb{R}^{I_{1}\times\ldots\times I_{N}\times
J_{1}\times\ldots\times J_{M}}$ and
$\mathcal{B}\in\mathbb{R}^{J_{1}\times\ldots\times J_{M}\times
K_{1}\times\ldots\times K_{L}}$ be tensors of orders $N+M$ and $M+L$,
respectively. Then

${\rm mat}(\mathcal{A}\ast_{M}\mathcal{B})={\rm mat}(\mathcal{A})\cdot{\rm
mat}(\mathcal{B}),$

where $\cdot$ denotes the usual matrix product.

Proof: The result follows by direct computations; see [9] for details. $\Box$

###### Definition 3

The following definitions can be found in, e.g., [9, 27, 31].

  1. 1.

The trace of a square tensor,
$\mathcal{A}\in\mathbb{R}^{I_{1}\times\ldots\times I_{N}\times
I_{1}\times\ldots\times I_{N}}$, is given by

${\rm
tr}(\mathcal{A})=\sum_{i_{1},\ldots,i_{N}}\mathcal{A}_{i_{1},\ldots,i_{N},i_{1},\ldots,i_{N}}.$

  2. 2.

The inner product of two tensors of the same size
$\mathcal{X},\mathcal{Y}\in\mathbb{R}^{I_{1}\times\ldots\times I_{N}\times
J_{1}\times\ldots\times J_{M}}$ is defined as

$\langle\mathcal{X},\mathcal{Y}\rangle=\sum_{i_{1},\ldots,i_{N},j_{1},\ldots,j_{M}}\mathcal{X}_{i_{1},\ldots,i_{N},j_{1},\ldots,j_{M}}\mathcal{Y}_{i_{1},\ldots,i_{N},j_{1},\ldots,j_{M}}.$

For square tensors, we have

$\langle\mathcal{X},\mathcal{Y}\rangle={\rm
tr}(\mathcal{X}^{T}\ast_{N}\mathcal{Y}).$

  3. 3.

The Frobenius norm of a tensor is given by

$\left|\left|\mathcal{X}\right|\right|_{F}=\sqrt{\langle\mathcal{X},\mathcal{X}\rangle}.$

If $\mathcal{X}$ is a square tensor, then

$\left|\left|\mathcal{X}\right|\right|_{F}=\sqrt{{\rm
tr}(\mathcal{X}^{T}\ast_{N}\mathcal{X})}.$

  4. 4.

For a positive integer $p$, we define the $p^{\rm th}$ power of a square
tensor $\mathcal{A}\in\mathbb{R}^{I_{1}\times\ldots\times I_{N}\times
I_{1}\times\ldots\times I_{N}}$ using the Einstein product recursively as

$\mathcal{A}^{p}=\mathcal{A}\ast_{N}\mathcal{A}^{p-1},$

where $\mathcal{A}^{0}=\mathcal{I}$ is the identity tensor.

###### Proposition 2

One has

$\left|\left|\mathcal{A}\right|\right|_{F}=\left|\left|{\rm
mat}(\mathcal{A})\right|\right|_{F}$ (1)

and, if $\mathcal{A}$ is a square tensor, then

$\left|\left|\mathcal{A}^{p}\right|\right|_{F}\leq\left|\left|\mathcal{A}\right|\right|_{F}^{p}.$
(2)

Proof: Let $\mathcal{A}\in\mathbb{R}^{I_{1}\times\ldots\times I_{N}\times
J_{1}\times\ldots\times J_{M}}$ be a tensor of order $N+M$. One has

$\left|\left|\mathcal{A}\right|\right|_{F}=\sqrt{\langle\mathcal{A},\mathcal{A}\rangle}=\sqrt{\sum_{i_{1},\ldots,i_{N},j_{1},\ldots,j_{M}}\mathcal{A}_{i_{1},\ldots,i_{N},j_{1},\ldots,j_{M}}^{2}}=\sqrt{\sum_{i=1}^{N}\sum_{j=1}^{M}A_{i,j}^{2}}=\left|\left|A\right|\right|_{F},$

with $A={\rm mat}(\mathcal{A})$ given in Definition 2. This shows (1). Assume
now that $M=N$. According to Definition 3(4.), by applying Proposition 1 $p-1$
times, one has ${\rm mat}(\mathcal{A}^{p})=({\rm mat}(\mathcal{A}))^{p}$, so
that $\left|\left|{\rm
mat}(\mathcal{A}^{p})\right|\right|_{F}\leq\left|\left|{\rm
mat}(\mathcal{A})\right|\right|_{F}^{p}$, where the inequality is due to the
submultiplicativity of the Frobenius matrix norm. Thus, thanks to (1), one has
inequality (2). This concludes the proof. $\Box$

###### Remark 4

The eigenvalue problem for $4^{\rm th}$-order tensors is discussed in [11].
The tensor format of this problem using the Einstein product is given by

$\mathcal{A}\ast_{2}X=\lambda X,$ (3)

where $\mathcal{A}\in\mathbb{R}^{I_{1}\times I_{2}\times J_{1}\times J_{2}}$,
$X\in\mathbb{R}^{J_{1}\times J_{2}}$ and $\lambda\in\mathbb{C}$. This
eigenvalue problem is equivalent to the matrix eigenvalue problem ${\rm
mat}(\mathcal{A})\cdot X=\lambda X$.

##  3 Tensor functions and centrality measures for multilayer networks

Centrality measures have been thoroughly studied for single-layer networks.
These measures include Katz centrality, subgraph centrality, and total
communicability with respect to a node; see, e.g., [5, 13, 17, 18, 19, 20]. In
this section, we introduce analogues of these measures for multilayer networks
using the adjacency tensor and the Einstein product. We limit our discussion
in this and the following sections to unweighted multilayer networks with only
one kind of connections in each layer. Thus, the adjacency tensor is a square
$4^{\rm th}$-order tensor. The extension to general weighted multilayer
networks, i.e., to multilayer networks with more than one kind of connections
in each layer, and with positive node-weights, is straightforward.

We recall that, for single-layer networks, a walk from node $v_{i_{1}}$ to
node $v_{i_{j}}$ is defined as a sequence of edges

$v_{i_{1}}\rightarrow v_{i_{2}},~{}v_{i_{2}}\rightarrow
v_{i_{3}},~{}\ldots,~{}v_{i_{j-1}}\rightarrow v_{i_{j}}$

that can be traversed to reach node $v_{i_{j}}$ from node $v_{i_{1}}$. The
length of the walk is $j-1$, the number of edges. Nodes and edges may be
repeated in a walk. A walk is said to be short if $j-1$ is fairly small.

In multilayer networks, there may be an edge from node $i$ in layer $\ell$ to
node $j$ in layer $k$. We denote this edge by $v_{i}^{\ell}\rightarrow
v_{j}^{k}$. A walk in a multilayer network from node $v_{i_{1}}^{\ell_{1}}$ to
node $v_{i_{j}}^{\ell_{j}}$ is defined as a sequence of edges

$v_{i_{1}}^{\ell_{1}}\rightarrow
v_{i_{2}}^{\ell_{2}},~{}v_{i_{2}}^{\ell_{2}}\rightarrow
v_{i_{3}}^{\ell_{3}},~{}\ldots,~{}v_{i_{j-1}}^{\ell_{j-1}}\rightarrow
v_{i_{j}}^{\ell_{j}}$

that can be traversed to reach node $v_{i_{j}}^{\ell_{j}}$ from node
$v_{i_{1}}^{\ell_{1}}$. The length of the walk is the number of edges, $j-1$,
and nodes and edges may be repeated in a walk. A closed multilayer walk is a
multilayer walk for which the starting and ending nodes are the same, i.e.,
$i_{1}=i_{j}$ and $\ell_{1}=\ell_{j}$. Estrada [16] has defined a walk in a
multiplex network in an analogous fashion.

The entries of the adjacency tensor $\mathcal{A}\in\mathbb{R}^{N\times L\times
N\times L}$ of an unweighted multilayer network tell us whether there is an
edge between any pairs of nodes (from the same or different layers). The
entries of the Einstein product of the adjacency tensor $\mathcal{A}$ with
itself,

$\mathcal{B}_{i,\ell,j,k}=\mathcal{A}\ast_{2}\mathcal{A}=\sum_{\begin{subarray}{c}1\leq
p\leq N\\\ 1\leq q\leq
L\end{subarray}}\mathcal{A}_{i,\ell,p,q}\mathcal{A}_{p,q,j,k},\qquad 1\leq
i,j\leq N,~{}~{}1\leq k,\ell\leq L,$

show the number of multilayer walks of length $2$ between pairs of nodes
$v_{i}^{\ell}$ and $v_{j}^{k}$. Similarly, let $p$ be a positive integer. Then
the entries of the tensor
$\mathcal{A}^{p}=\mathcal{A}\ast_{2}\mathcal{A}^{p-1}$ display the number of
multilayer walks of length $p$ between pairs of nodes $v_{i}^{\ell}$ and
$v_{j}^{k}$. In addition, the entry $\mathcal{A}^{p}_{i,\ell,i,\ell}$ provides
the number of closed multilayer walks of length $p$ that start at node
$v_{i}^{\ell}$. This suggests the introduction of the tensor function

$f(\mathcal{A})=\sum_{p=1}^{\infty}c_{p}\mathcal{A}^{p},$ (4)

where the coefficients $c_{p}$ generally are real and nonnegative, and are
chosen so that the series converges.

Let all entries of the matrix $E_{s,t}\in\mathbb{R}^{N\times L}$ vanish except
for the $(s,t)^{\rm th}$ entry, which is one. We refer to the entry

$f(\mathcal{A})_{i,\ell,j,k}=E_{i,\ell}\ast_{2}f(\mathcal{A})\ast_{2}E_{j,k}$
(5)

as the _communicability_ from node $v_{i}^{\ell}$ to node $v_{j}^{k}$; a
relatively large value indicates that it is easy to send information from node
$v_{i}^{\ell}$ to node $v_{j}^{k}$. Moreover, we call the entry

$f(\mathcal{A})_{i,\ell,i,\ell}=E_{i,\ell}\ast_{2}f(\mathcal{A})\ast_{2}E_{i,\ell}$
(6)

the _subgraph centrality_ of node $v_{i}^{\ell}$; a relatively large value
indicates that this node is important because much information may pass
through it. These notions of communicability and subgraph centrality for nodes
in multilayer graphs generalize the analogous definitions introduced and
explored by Estrada and his collaborators, as well as others, in [1, 13, 18,
19, 20] for single-layer networks. In these references $\mathcal{A}$ in (4) is
replaced by the adjacency matrix for the single-layer graph and, hence, $f$ is
a matrix function.

In many network applications short walks are more important than long walks,
because it is easier to transmit information via a few edges than via many
edges. This suggests that the coefficients $c_{p}$ should satisfy $0\leq
c_{p+1}\leq c_{p}$ for large $p$-values. One of the most commonly used matrix
functions for single-layer networks is the matrix exponential. We therefore
for multilayer networks introduce the tensor exponential

$\exp(\beta\mathcal{A})=\sum_{p=0}^{\infty}\frac{\beta^{p}\mathcal{A}^{p}}{p!}.$

It follows from Propositions 1 and 2 that this series converges for any fixed
$\beta$ in the interval $0\leq\beta<\infty$. The influence of the choice of
$\beta$ on the “$\delta$-effective diameter” on the coefficients $c_{p}$ of
the tensor function (4) is defined as the smallest integer $k$ such that

$\frac{\max_{\ell>k}|c_{\ell}|}{\max_{1\leq j\leq k}|c_{j}|}\leq\delta,$ (7)

for some small $\delta>0$. Its importance for the communicability in a network
is explored in [1, Definition 1]. Roughly, this diameter is $k$ if the tensor
function (4) can be approximated well by a polynomial of degree $k$. This
means that walks of length larger than $k$ do not significantly affect the
properties of the network. The definition in [1] differs slightly from (7) and
is for matrix functions, but carries over in a straightforward way to tensor
functions.

Since the first term $\mathcal{I}$ has no natural interpretation in the
context of network modeling, we will use the modified tensor exponential

$\exp_{0}(\mathcal{A}):=\exp(\mathcal{A})-\mathcal{I},$ (8)

where $\mathcal{I}$ denotes the identity tensor.

Resolvents of the adjacency matrix also are commonly used to determine
properties of nodes in a single-layer network; see, e.g., Estrada and Higham
[19] and Katz [26]. We define the modified tensor resolvent,

${\rm
res}_{0}(\mathcal{A},\alpha)=(\mathcal{I}-\alpha\mathcal{A})^{-1}-\mathcal{I}=\sum_{p=1}^{\infty}\alpha^{p}\mathcal{A}^{p},$
(9)

which is convergent for $0<\alpha<1/|\lambda_{\max}|$, where $\lambda_{\max}$
denotes eigenvalue of largest magnitude of $\mathcal{A}$. For many adjacency
tensors of interest, $\lambda_{\max}$ is real and positive. Conditions under
which this is the case are discussed by Qi and Luo [31]. For tensors of order
$4$, the eigenvalue $\lambda_{\max}$ can be computed as an eigenvalue of a
matrix using the relation of Remark 4. The choice of $\alpha$ affects the
$\delta$-effective diameter of the tensor function (9). This is discussed for
the matrix resolvent in [1].

We also define the _multilayer total communicability of node_ $v_{i}^{\ell}$
by

$E_{i,\ell}\ast_{2}f(\mathcal{A})\ast_{2}E,$ (10)

where all entries of the matrix $E\in\mathbb{R}^{N\times L}$ equal one, and
the _multilayer total communicability_ by

$E\ast_{2}f(\mathcal{A})\ast_{2}E.$ (11)

The latter definitions generalize analogous notions introduced for single-
layer networks by Benzi and Klymko [5] and Katz [26] for $f$ being the matrix
exponential or a matrix resolvent.

For small networks, we can evaluate the tensor functions discussed above by
applying the flattening operator ${\rm mat}$, its inverse ${\rm mat}^{-1}$,
Proposition 1, and using the following result.

###### Proposition 3

Let the tensor function $f$ be defined by (4) with a convergent power series.
Then

$f(\mathcal{A})={\rm mat}^{-1}(f({\rm mat}(\mathcal{A}))).$

Proof: By Proposition 1, ${\rm mat}(\mathcal{A}^{p})=({\rm
mat}(\mathcal{A}))^{p}$. Hence,

${\rm mat}(f(\mathcal{A}))={\rm
mat}(\sum_{p=1}^{\infty}c_{p}\mathcal{A}^{p})=\sum_{p=1}^{\infty}c_{p}{\rm
mat}(\mathcal{A}^{p})=\sum_{p=1}^{\infty}c_{p}({\rm mat}(\mathcal{A}))^{p}.$

Thus, one has

${\rm mat}(f(\mathcal{A}))=f({\rm mat}(\mathcal{A})).$

Applying the inverse operator ${\rm mat}^{-1}$ to both sides concludes the
proof. $~{}~{}~{}\Box$

The evaluation of tensor functions using the above proposition is feasible for
tensors that represent small to medium-sized multilayer networks. However, the
computations are very demanding for large-scale multilayer networks.
Approximations of tensor functions for the latter kind of networks can be
computed fairly inexpensively by applying Krylov subspace methods. This is
discussed in the following section.

##  4 Krylov subspace methods

Krylov subspace methods are well suited to approximate many matrix functions;
see, e.g., [2, 22] for illustrations. They also have been applied successfully
to the approximation of tensor functions and the solution of tensor systems of
equations; see [3, 4, 14, 23, 30] and references therein. It is therefore
natural to seek to approximate the tensor functions mentioned in the previous
section by Krylov subspace methods. We first discuss the application of the
global tensor Arnoldi process to the approximation of multilayer centrality
measures, and subsequently consider the tensor block Arnoldi process.

###  4.1 A global tensor Arnoldi process based on the Einstein product

The global matrix Arnoldi process is a Krylov subspace method that was
introduced by Jbilou et al. [24, 25] for the reduction of a large matrix to a
small one. A global tensor Arnoldi process for the reduction of a large tensor
to a small matrix using the Einstein product is described by El Guide et al.
[14]. The application of $m$ steps of the latter process to the tensor
$\mathcal{A}\in\mathbb{R}^{N\times L\times N\times L}$ with initial matrix
$V\in\mathbb{R}^{N\times L}$ determines, when no breakdown occurs, an
orthonormal basis for the tensor Krylov subspace

$\mathcal{K}_{m+1}(\mathcal{A},V)={\rm
span}\\{V,\mathcal{A}\ast_{2}V,\ldots,\mathcal{A}^{m}\ast_{2}V\\}:=\left\\{\sum_{i=0}^{m}\omega_{i}\mathcal{A}^{i}\ast_{2}V,~{}\omega_{i}\in\mathbb{R}\right\\}.$
(12)

Algorithm 1 Global tensor Arnoldi process

  * 1:

Input: Adjacency tensor $\mathcal{A}\in\mathbb{R}^{{N\times L\times N\times
L}}$, initial block vector $V\in\mathbb{R}^{N\times L}$, and number of steps
$m$.

  * 2:

Output: Orthonormal basis $\mathcal{V}_{m+1}=\\{V_{1},V_{2},\ldots,V_{m+1}\\}$
for the tensor Krylov subspace (12) and nontrivial entries of the Hessenberg
matrix $H_{m+1,m}=[h_{ij}]\in\mathbb{R}^{(m+1)\times m}$.

  * 3:

$V_{1}=V/\|V\|_{F}$

  * 4:

for $j=1,\ldots,m$ do

  * 5:

$W=\mathcal{A}\ast_{2}V_{j}$

  * 6:

for $i=1,\ldots,j$ do

  * 7:

$h_{i,j}=(V_{i},W)$

  * 8:

$W=W-h_{ij}V_{j}$

  * 9:

end for

  * 10:

$h_{j+1,j}=\left|\left|W\right|\right|_{F}$

  * 11:

if $h_{j+1,j}=0$, then stop

  * 12:

else $V_{j+1}=W/h_{j+1,j}$

  * 13:

end if

  * 14:

end for

The computations are described by Algorithm 1. The algorithm is said to _break
down_ at step $j$ if $h_{i+1,i}>0$ for $1\leq i<j$ and $h_{j+1,j}=0$. In the
absence of breakdown, the algorithm determines the tensor
$\mathcal{V}_{m+1}=[V_{1},V_{2},\ldots,V_{m+1}]\in\mathbb{R}^{N\times
L\times(m+1)}$ with orthonormal block columns, i.e.,

$(V_{i},V_{j}):={\rm tr}(V_{i}^{T}V_{j})=\left\\{\begin{array}[]{cc}1,&i=j,\\\
0,&i\neq j,\end{array}\right.$

that span the tensor Krylov subspace $\mathcal{K}_{m+1}(\mathcal{A},V)$. In
line 5 of Algorithm 1, we have

$\mathcal{A}\ast_{2}\mathcal{V}_{j}=[\mathcal{A}\ast_{2}V_{1},\mathcal{A}\ast_{2}V_{2},\ldots,\mathcal{A}\ast_{2}V_{j}]\in\mathbb{R}^{N\times
L\times j}.$

It follows from the recursion relation of Algorithm 1 that

$\mathcal{A}\ast_{2}\mathcal{V}_{m}=\mathcal{V}_{m+1}\ast_{1}H_{m+1,m},$ (13)

where $H_{m+1,m}=[h_{ij}]\in\mathbb{R}^{(m+1)\times m}$ is an upper Hessenberg
matrix made up of the coefficients $h_{ij}$ generated in lines 7 and 10 of
Algorithm 1; all entries below the subdiagonal of $H_{m+1,m}$ vanish.

Let the matrix $H_{m}\in\mathbb{R}^{m\times m}$ be obtained by deleting the
last row of $H_{m+1,m}$. Then

$\mathcal{V}_{m}^{T}\ast_{2}\mathcal{A}\ast_{2}\mathcal{V}_{m}=H_{m},$

where $\mathcal{V}_{m}=[V_{1},V_{2},\ldots,V_{m}]\in\mathbb{R}^{N\times
L\times m}$ and
$\mathcal{V}_{m}^{T}=[V_{1},V_{2},\ldots,V_{m}]^{T}\in\mathbb{R}^{m\times
N\times L}$. Hence, $H_{m}$ is the orthogonal projection of $\mathcal{A}$ onto
the subspace $\mathcal{K}_{m}(\mathcal{A},V)$ with respect the basis
$\mathcal{V}_{m}$. This suggest to use the approximation

$\mathcal{V}_{m}\ast_{1}f(H_{m})\ast_{1}E_{1}\|V\|_{F}$ (14)

of $f(\mathcal{A})\ast_{2}V$, where $\|V\|_{F}=\sqrt{{\rm tr}(V^{T}V)}$,
analogously to the approach used when $\mathcal{A}$ is a square matrix and $V$
is a vector; see [2, 14, 23].

This approach to approximate $f(\mathcal{A})\ast_{2}V$ works well when $V=E$
and can be applied to determine accurate approximations of the multilayer
total communicability (11) and the multilayer total communicability of node
$v_{i}^{\ell}$ defined by (10). The former is approximated by

$E\ast_{2}\mathcal{V}_{m}\ast_{1}f(H_{m})\ast_{1}E_{1}\|V\|_{F}$

and the latter by

$E_{i,\ell}\ast_{2}\mathcal{V}_{m}\ast_{1}f(H_{m})\ast_{1}E_{1}\|V\|_{F}.$
(15)

In particular, the evaluation of (15) does not require any arithmetic work
when the expression (14) is available. This makes the evaluation of the
multilayer total communicability of all nodes $v_{i}^{\ell}$, $1\leq i\leq N$
and $1\leq\ell\leq L$, inexpensive when the expression (14) is known. We use
this fact when determining nodes for which this measure is large in Section 5.

However, Algorithm 1 often suffers from breakdown when seeking to approximate
an expression of the form $f(\mathcal{A})\ast_{2}V$ when the matrix $V$ is
sparse, i.e., when $V$ has many vanishing entries. This is the case when
seeking to approximate the subgraph centrality (6) by

$E_{i,\ell}\ast_{2}\mathcal{V}_{m}\ast_{1}f(H_{m})\ast_{1}E_{1}\|E_{i,\ell}\|_{F},$

or the communicability (5) between the node $v_{i}^{\ell}$ and node
$v_{j}^{k}$ by

$E_{i,\ell}\ast_{2}\mathcal{V}_{m}\ast_{1}f(H_{m})\ast_{1}E_{1}\|E_{j,k}\|_{F}.$

When computing these approximations, the initial block matrix is
$V=E_{i,\ell}$, which is very sparse. Since the tensor $\mathcal{A}$ typically
also is sparse, this often results in that the scalar $h_{j+1,j}$ in line 10
of Algorithm 1 vanishes for some $1\leq j\leq m$. The computations with the
algorithm then cannot be continued, and the available expression at breakdown,

$\mathcal{V}_{j}\ast_{1}f(H_{j})\ast_{1}E_{1}\|V\|_{F},$

might not furnish an approximation of desired accuracy. Moreover, even when
all the multilayer subgraph centralities can be computed to determine the node
with the largest subgraph centrality, this is quite expensive for large
multilayer networks. We describe in the following subsection how these
difficulties can be reduced by replacing the initial block vector $V$ in
Algorithm 1 by a tensor.

###  4.2 A block Arnoldi process based on the Einstein product

We describe a block Arnoldi process that uses the Einstein product. It differs
from Algorithm 1 in that the initial block vector $V$ is extended to a tensor.
The application of $m$ steps of the block Arnoldi process to $\mathcal{A}$
with initial tensor $\mathcal{V}$ determines, in the absence of breakdown, an
orthonormal basis for the block tensor Krylov subspace

$\mathcal{K}_{m+1}^{\rm block}(\mathcal{A},\mathcal{V})={\rm
span}\\{\mathcal{V},\mathcal{A}\ast_{2}\mathcal{V},\ldots,\mathcal{A}^{m}\ast_{2}\mathcal{V}\\}=\left\\{\sum_{i=0}^{m}\mathcal{A}^{i}\ast_{2}\mathcal{V}\ast_{1}\Omega_{i},~{}\Omega_{i}\in\mathbb{R}^{R\times
R}\right\\},$ (16)

where $\mathcal{V}\in\mathbb{R}^{N\times L\times R}$. We will refer to the
integer $R$ as the block size. The block Arnoldi process of this subsection
has the advantage of typically reducing the number of times we need to access
the adjacency tensor when compared with Algorithm 1. Moreover, choosing
suitable auxiliary columns in the initial tensor $\mathcal{V}$, the
occurrences of breakdowns can be reduced in comparison with Algorithm 1. The
block Arnoldi process is summarized by Algorithm 2.

Algorithm 2 Block tensor Arnoldi process

  * 1:

Input: Adjacency tensor $A\in\mathbb{R}^{{N\times L\times N\times L}}$,
initial tensor $\mathcal{V}\in\mathbb{R}^{N\times L\times R}$, and number of
steps $m$.

  * 2:

Output: Orthonormal basis
$\mathbb{V}_{m+1}=\\{\mathcal{V}_{1},\mathcal{V}_{2},\ldots,\mathcal{V}_{m+1}\\}$
for the block Krylov subspace (16), and nontrivial entries
$H_{i,j}\in\mathbb{R}^{R\times R}$ of the block Hessenberg matrix
$\mathbb{H}_{m+1,m}=[H_{i,j}]\in\mathbb{R}^{R(m+1)\times Rm}$.

  * 3:

Compute the QR factorization of $\mathcal{V}=\mathcal{Q}\ast_{1}R$, where the
tensor $\mathcal{Q}\in\mathbb{R}^{N\times L\times R}$ satisfies
$\mathcal{Q}^{T}\ast_{1}\mathcal{Q}=\mathcal{I}$ and the matrix
$R\in\mathbb{R}^{R\times R}$ is upper triangular. Set
$\mathcal{V}_{1}=\mathcal{Q}$ and $H_{1,0}=R$.

  * 4:

for $j=1,\ldots,m$ do

  * 5:

$\mathcal{W}=\mathcal{A}\ast_{2}\mathcal{V}_{j}$

  * 6:

for $i=1,\ldots,j$ do

  * 7:

$H_{i,j}=\mathcal{V}_{i}^{T}\ast_{2}\mathcal{W},$

  * 8:

$\mathcal{W}=\mathcal{W}-\mathcal{V}_{i}\ast_{1}H_{i,j}$

  * 9:

end for

  * 10:

Compute the QR factorization $\mathcal{W}=\mathcal{Q}\ast_{1}R$, where
$\mathcal{Q}\in\mathbb{R}^{N\times L\times R}$ satisfies
$\mathcal{Q}^{T}\ast_{1}\mathcal{Q}=\mathcal{I}$ and the matrix
$R\in\mathbb{R}^{R\times R}$ is upper triangular. Set
$\mathcal{V}_{j+1}=\mathcal{Q}$ and $H_{j+1,j}=R$.

  * 11:

end for

Algorithm 2 determines an orthonormal basis

$\mathbb{V}_{m+1}=\left[\mathcal{V}_{1},\ldots,\mathcal{V}_{m+1}\right]\in\mathbb{R}^{N\times
L\times R(m+1)}$

for the block Krylov subspace (16) and the upper block Hessenberg matrix
$\mathbb{H}_{m+1,m}=[H_{i,j}]\in\mathbb{R}^{R(m+1)\times Rm}$. We have the
following result:

###### Proposition 4

Suppose that $m$ steps of Algorithm 2 have been carried out. Then

$\mathcal{A}\ast_{2}\mathbb{V}_{m}=\mathbb{V}_{m+1}\ast_{1}\mathbb{H}_{m+1,m}$
(17)

and

$\mathcal{A}\ast_{2}\mathbb{V}_{m}=\mathbb{V}_{m}\ast_{1}\mathbb{H}_{m}+\mathcal{V}_{m+1}\ast_{1}H_{m+1,m}\ast_{1}\mathbb{E}_{m}^{T},$
(18)

where
$\mathbb{V}_{m}=[\mathcal{V}_{1},\ldots,\mathcal{V}_{m}]\in\mathbb{R}^{N\times
L\times Rm}$ is made up of the first $m$ tensors of $\mathbb{V}_{m+1}$,
$\mathbb{E}_{m}^{T}=[0,0,\ldots,I_{R}]\in\mathbb{R}^{R\times Rm}$ and

$\mathbb{H}_{m+1,m}=\begin{bmatrix}&&\mathbb{H}_{m}&\\\
0&\cdots&0&H_{m+1,m}\end{bmatrix}.$

Moreover,

$\mathcal{A}^{p}\ast_{2}\mathcal{V}=\mathbb{V}_{m}\ast_{1}\mathbb{H}_{m}^{p}\ast_{1}E_{1}\chi_{0},$
(19)

where $\chi_{0}$ is such that
$\mathcal{V}=\mathbb{V}_{m}\ast_{1}E_{1}\chi_{0}$ holds for all
$p=0,1,\dots,m-1$.

Proof: For $1\leq j\leq m$, one has

$\displaystyle[\mathbb{V}_{m+1}\ast_{1}\mathbb{H}_{m+1,m}]_{:,:,1+R(j-1):Rj}$
$\displaystyle=\sum_{i=1}^{j+1}(\mathbb{V}_{m+1})_{:,:,1+R(i-1):Ri}\ast_{1}(\mathbb{H}_{m+1,m})_{1+R(i-1):Ri,1+R(j-1):Rj}$
$\displaystyle=[\mathcal{A}\ast_{2}\mathbb{V}_{m}]_{:,:,1+R(j-1):Rj}.$

This leads to Equation (17); Equation (18) can be shown similarly. The last
claim can be proved by induction. In fact, if $p=0$, then one has
$\mathcal{A}^{0}\ast_{2}\mathcal{V}=\mathcal{V}=\mathbb{V}_{m}\ast_{1}E_{1}\chi_{0}$
and, by assuming that (19) holds for $p=0,1,\ldots,m-2$, one obtains

$\mathcal{A}^{p+1}\ast_{2}\mathcal{V}=\mathcal{A}\ast_{2}\mathcal{A}^{p}\ast_{2}\mathcal{V}=\mathcal{A}\ast_{2}\mathbb{V}_{m}\ast_{1}\mathbb{H}_{m}^{p}\ast_{1}E_{1}\chi_{0},$

so that, using Equation (18),

$\displaystyle\mathcal{A}^{p+1}\ast_{2}\mathcal{V}$
$\displaystyle=(\mathbb{V}_{m}\ast_{1}\mathbb{H}_{m}+\mathcal{V}_{m+1}\ast_{1}H_{m+1,m}\ast_{1}\mathbb{E}_{m}^{T})\ast_{1}\mathbb{H}_{m}^{p}\ast_{1}E_{1}\chi_{0}$
$\displaystyle=\mathbb{V}_{m}\ast_{1}\mathbb{H}_{m}\ast_{1}\mathbb{H}_{m}^{p}\ast_{1}E_{1}\chi_{0}+\mathcal{V}_{m+1}\ast_{1}H_{m+1,m}\ast_{1}\mathbb{E}_{m}^{T}\ast_{1}\mathbb{H}_{m}^{p}\ast_{1}E_{1}\chi_{0},$

where the second term vanishes due to the fact that $\mathbb{H}_{m}$ is a
Hessenberg matrix. This concludes the proof. $~{}~{}~{}\Box$

Thanks to Equation (19), $f(\mathcal{A})\ast_{2}\mathcal{V}$ can be
approximated by
$\mathbb{V}_{m}\ast_{2}f(\mathbb{H}_{m})\ast_{1}E_{1}\chi_{0}$, where
$\chi_{0}$ is such that $\mathcal{V}=\mathbb{V}_{m}\ast_{1}E_{1}\chi_{0}$. The
advantage of this approach is that we can compute the multilayer subgraph
centrality and the resolvent-based subgraph centrality of $R$ nodes at once,
and we also can determine approximations of the multilayer communicabilities
of these $R$ nodes essentially for free. This is because we are approximating
the quantity $f(\mathcal{A})\ast_{2}\mathcal{V}$, and then can evaluate
$\mathcal{V}\ast_{2}f(\mathcal{A})\ast_{2}\mathcal{V}$ inexpensively.
Moreover, we can always correct the numerical stability issue related to the
sparsity of the adjacency tensor and the initial block by adding a dense
matrix-column in the initial block. This technique has been discussed for
matrix functions case in [21]. Notice that if we include $E$, i.e., the
matrix-column of all ones, in the initial block, then Algorithm 2 will produce
the same quantities of interest as one would compute with Algorithm 1 when
applied with an initial matrix $V=E$.

##  5 Computed examples

This section presents some examples to illustrate the performance of the
methods discussed above. The computations were carried out using MATLAB
R2015b. For the examples in Sections 5.2 and 5.3, we choose the minimum number
of steps, $m$, with the Krylov subspace method needed to obtain the same
ranking of the first $10$ nodes as the ranking obtained when evaluating the
exact tensor function. We will see that the number of steps required is quite
small. Due to the size of the networks in the examples in Sections 5.4 and
5.5, it is expensive to evaluate the exact tensor function. In these examples,
we therefore increase the number of steps, $m$, of the Krylov subspace until
the ranking does not change.

$1$$2$$3$$4$$5$$1$$2$$3$$4$$5$ Figure 1: Example 1: Layers are presented from
left to right in the order $L=1$ and $L=2$. The edges connecting nodes from
same layer are marked in black. The edges connecting nodes from different
layers are marked in red.

###  5.1 Example 1: A small synthetic multilayer network

Regard a small synthetic unweighted and undirected multilayer network with
$L=2$ layers and $N=5$ nodes in each layer. It consists of $10$ edges as shown
by the graph of Figure 1. Let $\mathcal{A}$ be the adjacency tensor. We
compute the different multilayer centrality measures discussed using the
Einstein product. The results are summarized in Table 1 for $\beta=1$ and
$\alpha=0.5/\lambda_{max}$. We can see that node $v_{2}^{2}$ has the largest
centrality measures and, therefore, is the most node, i.e., the most important
node. This node is followed by node $v_{1}^{1}$. Node $v_{4}^{2}$ is least
important having only one connection. All centrality measures perform as can
be expected.

Table 1: Multilayer total communicability of the nodes $v_{i}^{\ell}$ (MTC), multilayer Katz centrality (MKC), multilayer subgraph centrality with the modified tensor exponential (MSC${}_{\exp_{0}}$), multilayer subgraph centrality with the modified tensor resolvent (MSC${}_{\rm res_{0}}$) for the nodes of Example 1. | $MTC(i,\ell)$ |  | $MKC(i,\ell)$ |
---|---|---|---|---
$v_{i}^{\ell}$ | $\ell=1$ | $\ell=2$ | $\ell=1$ | $\ell=2$
$i=1$ | $12.2520$ | $7.7379$ | $2.1131$ | $1.7044$
$i=2$ | $9.6537$ | $17.2450$ | $1.8145$ | $2.5454$
$i=3$ | $8.9474$ | $11.7351$ | $1.7645$ | $1.9484$
$i=4$ | $10.8250$ | $4.2550$ | $1.8876$ | $1.3470$
$i=5$ | $10.6175$ | $10.6175$ | $1.8774$ | $1.8774$
| $MSC_{\exp_{0}}(i,\ell)$ |  | $MSC_{\rm res_{0}}(i,\ell)$ |
$v_{i}^{\ell}$ | $\ell=1$ | $\ell=2$ | $\ell=1$ | $\ell=2$
$i=1$ | $3.1001$ | $2.2834$ | $1.1507$ | $1.0952$
$i=2$ | $2.3582$ | $4.1313$ | $1.0987$ | $1.2165$
$i=3$ | $2.3946$ | $2.4698$ | $1.1003$ | $1.1039$
$i=4$ | $2.4174$ | $1.5922$ | $1.1016$ | $1.0454$
$i=5$ | $2.5001$ | $2.5001$ | $1.1051$ | $1.1051$

###  5.2 Example 2: A synthetic multilayer network

We consider an example of a weighted multilayer network with $L=32$ layers,
$N=20$ nodes in each layer, and $674$ edges. The network can be downloaded
from https://github.com/wjj0301/Multiplex-Networks. We compute approximations
of the multilayer total communicability of the nodes $v_{i}^{\ell}$ and
multilayer Katz centrality using Algorithm 1, as well as the exact multilayer
total communicability of the nodes $v_{i}^{\ell}$ and multilayer Katz
centrality by flattening the adjacency tensor and using the MATLAB function
expm and the MATLAB backslash operator. The errors in the computed
approximations is determined by the vector infinity norm, that is the exact
multilayer total communicability of the nodes, multilayer Katz centrality, and
their approximations are stored in vectors and the infinity norm is applied to
measure the distance between the vectors with the exact entries and the
vectors with the corresponding entries determined by Krylov subspace methods.
Figure 2 displays the errors as a function of the Krylov subspace dimension
$m$. Table 2 lists the top $10$ central nodes according to the multilayer
total communicability with respect to a node and multilayer Katz centrality.
We notice that the multilayer Katz centrality is approximated accurately when
the same number of steps, $m$, of Algorithm 1 are carried out. The multilayer
total communicability of the nodes is approximated less accurately with the
same number of steps, but still gives the same ranking of nodes as when the
exact communicability is used.

Table 2: Top $10$ central nodes determined by Algorithm 1 with $m=6$ using the (approximate) multilayer total communicability of the nodes (MTC) and (approximate) multilayer Katz centrality (MKC) so obtained for the weighted multilayer network in Example 2, with $\alpha=0.4/\lambda_{\max}$ and $\beta=0.4$. These computed values are compared to the exact ones obtained by flattening and using matrix functions. $\\{i,\ell\\}$ | $MTC\\{i,\ell\\}$ | $MTC_{\rm exact}\\{i,\ell\\}$ | $\\{i,\ell\\}$ | $MKC\\{i,\ell\\}$ | $MKC_{\rm exact}\\{i,\ell\\}$
---|---|---|---|---|---
$\\{18,24\\}$ | $99.4969$ | $97.1435$ | $\\{18,24\\}$ | $3.6474$ | $3.6407$
$\\{17,26\\}$ | $70.0610$ | $72.6593$ | $\\{19,19\\}$ | $3.0148$ | $3.0191$
$\\{13,19\\}$ | $64.5009$ | $67.7357$ | $\\{14,23\\}$ | $2.8966$ | $2.8956$
$\\{8,26\\}$ | $62.8084$ | $67.1222$ | $\\{13,19\\}$ | $2.6318$ | $2.6407$
$\\{19,19\\}$ | $62.1558$ | $63.7426$ | $\\{17,26\\}$ | $2.6242$ | $2.6313$
$\\{6,24\\}$ | $59.3319$ | $60.2469$ | $\\{6,29\\}$ | $2.5999$ | $2.5990$
$\\{19,4\\}$ | $58.9334$ | $59.4300$ | $\\{2,29\\}$ | $2.5497$ | $2.5477$
$\\{14,24\\}$ | $56.5841$ | $56.4157$ | $\\{1,32\\}$ | $2.5371$ | $2.5425$
$\\{1,32\\}$ | $52.8265$ | $54.7537$ | $\\{1,3\\}$ | $2.5251$ | $2.5240$
$\\{2,24\\}$ | $51.6585$ | $54.3109$ | $\\{14,24\\}$ | $2.5150$ | $2.5146$

We also set the weight of all edges to one to obtain an unweighted multilayer
network and compute the multilayer total communicability of the nodes and
multilayer Katz centrality. We notice that multilayer total communicability of
the nodes is approximated accurately when $m=6$ steps of Algorithm 1 are
carried out also for larger values of $\beta$. However, the multilayer Katz
centrality is not approximated accurately for $m=6$ steps when
$\alpha>0.66/\lambda_{\max}$. Figure 2 depicts the influence of the values of
$\alpha$ and $\beta$ on the dimension of Krylov subspace needed to obtain
accurate approximations for both weighted and unweighted graphs.

Finally, we turn to the computation of multilayer subgraph centralities of
some nodes. We apply Algorithm 2 to our adjacency tensor with the first block
$\mathcal{V}\in\mathbb{R}^{N\times L\times R}$ with $R$ frontal slices of the
form $E_{i,\ell}$. We let $R=10$ and choose the indices $\\{i,\ell\\}$
randomly. The results obtained are summarized in Table 3, which also shows
communicabilities that we obtain for free since we use a block method.

(a)

(b)

(c)

(d)

Figure 2: Infinity norm error for multilayer total communicability of the nodes and multilayer Katz centrality as functions of the Krylov subspace dimension $m$ for Example $2$ with $\alpha\in\\{\frac{0.2}{\lambda_{\max}},\frac{0.3}{\lambda_{\max}},\frac{0.5}{\lambda_{\max}},\frac{0.7}{\lambda_{\max}}\\}$, and $\beta\in\\{0.3,0.5,1,1.5\\}$ ($\\{$blue, red, yellow, green$\\}$). $(a)$ multilayer Katz centrality for weighted edges, $(b)$ multilayer total communicability of the nodes for weighted edges, $(c)$ multilayer Katz centrality for unweighted edges, and $(d)$ multilayer total communicability of the nodes for unweighted edges. Table 3: Multilayer subgraph centrality obtained with the modified tensor exponential (MSC${}_{\exp_{0}}$) and multilayer subgraph centrality with the modified tensor resolvent (MSC${}_{\rm res_{0}}$) for some nodes for the unweighted multilayer network in Example 2, as well as free multilayer communicabilities of the nodes (MC) with $\alpha=0.5/\lambda_{\max}$, $\beta=1$, and $m=5$. $\\{i,\ell\\}$ | $MSC_{\exp_{0}}\\{i,\ell\\}$ | $MSC_{\rm res_{0}}\\{i,\ell\\}$ | $\\{i,\ell,j,k\\}$ | $MC\\{i,\ell,j,k\\}$
---|---|---|---|---
$\\{4,1\\}$ | $1.1683$ | $1.0358$ |  |
$\\{5,1\\}$ | $1.2116$ | $1.0488$ |  |
$\\{12,1\\}$ | $1.5878$ | $1.1321$ |  |
$\\{1,3\\}$ | $1.1766$ | $1.0395$ | $\\{5,1,4,1\\}$ | $1.0505$
$\\{19,4\\}$ | $1.0084$ | $1.0038$ | $\\{12,1,4,1\\}$ | $0.1768$
$\\{2,11\\}$ | $1.0592$ | $1.0193$ | $\\{12,1,5,1\\}$ | $0.5531$
$\\{13,18\\}$ | $1.1681$ | $1.0354$ | $\\{13,18,19,18\\}$ | $0.0502$
$\\{19,18\\}$ | $1.0097$ | $1.0048$ |  |
$\\{14,22\\}$ | $1.2537$ | $1.0611$ |  |
$\\{8,25\\}$ | $1.5431$ | $1.1176$ |  |

###  5.3 Example 3: Multiplex network (Scotland Yard transportation data)

This example considers the Scotland Yard transportation network created by the
authors of [7], which is multiplex network. A multiplex network is a special
case of a multilayer network. The network can be downloaded from [6] as a
weighted or unweighted multiplex network. It consists of $3324$ edges and
$N=199$ nodes that represent public transport stops in the city of London. The
network has $L=4$ layers that represent different modes of transportation:
Boat, underground, bus, and taxi. The weights are determined so that the edges
in the layer that represent travel by taxi all have weight one. A taxi ride is
defined as a trip by taxi between two adjacent nodes in the taxi layer; a taxi
ride along $k$ edges is considered $k$ taxi rides. The weights of edges in the
boat, underground, and bus layers are chosen to be equal to the minimal number
of taxi rides required to travel between the same nodes. We compute the
multilayer Katz centrality and multilayer total communicability of the nodes
by applying Algorithm 1 to the adjacency tensor of the given multiplex
network, as well as to the adjacency tensor of the associated unweighted
network; in the latter adjacency tensor all edges have weight one. We compute
the exact multilayer Katz centrality and multilayer total communicability of
the nodes in the same way as in the previous example and evaluate the infinity
norm error as well. Figure 3 displays the errors as a function of the Krylov
subspace dimension. This illustrates the accuracy of Algorithm 1 when applied
to multiplex networks.

(a)

(b)

Figure 3: Infinity norm error for multilayer total communicability of the
nodes (blue) and multilayer Katz centrality (orange) as functions of the
Krylov subspace dimension $m$ for Example $3$. $(a)$ when edges are unweighted
and $\alpha=0.4/\lambda_{\max}$ and $\beta=0.4$, $(b)$ when edges are weighted
and $\alpha=0.3/\lambda_{\max}$ and $\beta=0.1$.

###  5.4 Example 4: Multiplex network (European airlines data set)

The European airlines data set consists of $N=450$ nodes that represent
European airports and has $L=37$ layers that represent different airlines
operating in Europe. There are $3588$ edges, which represent available routs.
This network can be represented by a fourth-order adjacency tensor
$\mathcal{A}\in\mathbb{R}^{N\times L\times N\times L}$ such that
$\mathcal{A}(i,\ell,j,\ell)=1$ if there is a flight connecting airports $i$
and $j$ with airline $\ell$. Moreover, $\mathcal{A}(i,\ell,i,k)=1$ for every
$1\leq\ell,k\leq L$ to reflect the effort required to change airlines for
connecting flights. The network can be downloaded from [6]. Similarly as
Taylor et al. [32], we only include $N=417$ nodes from the largest connected
component of the network. We compute the multilayer total communicability of
the nodes and the multilayer Katz centrality using Algorithm 1 to approximate
the tensor exponential and the tensor resolvent functions. Table 4 lists the
top $10$ central nodes. We obtained similar ranking as reported in [7], where
the authors applied Krylov subspace methods to the supra-adjacency matrix of
the network in order to compute matrix function-based centrality measures such
as the Katz centrality.

Table 4: Top $10$ central nodes according to multilayer total communicability of the nodes (MTC) and multilayer Katz centrality (MKC) for the European airlines network in Example 4, with $\alpha=0.5/\lambda_{\max}$ and $\beta=0.2$. The MTC and MKC values are computed with Algorithm 1 with $m=20$.  $\\{i,\ell\\}$ | $MTC\\{i,\ell\\}$ | $\\{i,\ell\\}$ | $MKC\\{i,\ell\\}$
---|---|---|---
$\\{\text{Stansted, Ryanair}\\}$ | $8.2164\cdot 10^{3}$ | $\\{\text{Stansted, Ryanair}\\}$ | $4.4228$
$\\{\text{Munich, Lufthansa}\\}$ | $7.5209\cdot 10^{3}$ | $\\{\text{Munich, Lufthansa}\\}$ | $4.0937$
$\\{\text{Frankfurt, Lufthansa}\\}$ | $7.4541\cdot 10^{3}$ | $\\{\text{Frankfurt, Lufthansa}\\}$ | $4.0650$
$\\{\text{Dublin, Ryanair}\\}$ | $6.8942\cdot 10^{3}$ | $\\{\text{Ataturk, Turkish}\\}$ | $4.0486$
$\\{\text{Gatwick, EasyJet}\\}$ | $6.5714\cdot 10^{3}$ | $\\{\text{Gatwick, EasyJet}\\}$ | $3.7925$
$\\{\text{Ataturk, Turkish}\\}$ | $6.4399\cdot 10^{3}$ | $\\{\text{Dublin, Ryanair}\\}$ | $3.6479$
$\\{\text{Amsterdam, KLM}\\}$ | $6.0994\cdot 10^{3}$ | $\\{\text{Vienna, Austrian}\\}$ | $3.5939$
$\\{\text{Vienna, Austrian}\\}$ | $5.8057\cdot 10^{3}$ | $\\{\text{Amsterdam, KLM}\\}$ | $3.5661$
$\\{\text{Caravaggio, Ryanair}\\}$ | $5.7104\cdot 10^{3}$ | $\\{\text{Caravaggio, Ryanair}\\}$ | $3.3244$
$\\{\text{Adolfo, Ryanair}\\}$ | $5.5765\cdot 10^{3}$ | $\\{\text{Charles de Gaulle, Air France}\\}$ | $3.2444$

###  5.5 Example $5$: Wikispeedia network

This dataset contains human navigation paths in Wikipedia, collected through
the human-computation game Wikispeedia. Wikispeedia users are asked to
navigate from a given source to a given target article, being allowed to click
on links only. Nodes are articles of the English Wikipedia and edges represent
clicks. The data is provided by the Stanford Network Analysis Project
http://snap.stanford.edu/index.html. It contains $4604$ articles and $119882$
links. We classified the articles into $16$ different subjects (Countries,
Science, Geography, …) and then built a multilayer network with $N=4604$ nodes
and $L=16$ layers. Each layer contains the edges connecting the nodes that are
considered to be classified in this layer. We have used the files that contain
node identifiers and all edges from https://github.com/franloza/Wikispeedia-
Network. We compute the multilayer Katz centrality and multilayer total
communicability of the nodes by applying Algorithm 1 to the adjacency tensor
of the given multilayer network. Table 5 lists the top $10$ central nodes for
this network.

Table 5: Top $10$ central nodes according to the multilayer total communicability (MTC) and the multilayer Katz centrality (MKC) for the Wikispeedia network in Example $5$, with $\beta=0.2$, $\alpha=0.5/\lambda_{\max}$. The MTC and MKC values are approximations determined by Algorithm 1 with $m=20$. $\\{i,\ell\\}$ | $MTC\\{i,\ell\\}$
---|---
$\\{\text{Fauna of Australia, Science}\\}$ | $2.3711\cdot 10^{27}$
$\\{\text{Africa, Countries}\\}$ | $2.1895\cdot 10^{27}$
$\\{\text{Periodic table, Science}\\}$ | $1.7854\cdot 10^{27}$
$\\{\text{List of elements by name, Science}\\}$ | $1.7241\cdot 10^{27}$
$\\{\text{Periodic table (large version), Science}\\}$ | $1.7162\cdot 10^{27}$
$\\{\text{Bird, Science}\\}$ | $1.5887\cdot 10^{27}$
$\\{\text{Bird, Animals}\\}$ | $1.5887\cdot 10^{27}$
$\\{\text{President of the United States, Politics}\\}$ | $1.5476\cdot 10^{27}$
$\\{\text{Star, Science}\\}$ | $1.5465\cdot 10^{27}$
$\\{\text{Star, Space}\\}$ | $1.5465\cdot 10^{27}$
$\\{i,\ell\\}$ | $MKC\\{i,\ell\\}$
$\\{\text{Lebanon, Geography}\\}$ | $4.7019$
$\\{\text{Armenia, Countries}\\}$ | $4.5967$
$\\{\text{Armenia, Geography}\\}$ | $4.5967$
$\\{\text{Georgia, Countries}\\}$ | $4.5712$
$\\{\text{Georgia, Geography}\\}$ | $4.5712$
$\\{\text{Turkey, Countries}\\}$ | $4.4296$
$\\{\text{Turkey, Geography}\\}$ | $4.4296$
$\\{\text{Djibouti, Countries}\\}$ | $4.2680$
$\\{\text{Djibouti, Geography}\\}$ | $4.2680$
$\\{\text{Mozambique, Countries}\\}$ | $4.1781$

##  6 Conclusion

This paper investigates centrality measures for multilayer networks by
introducing the exponential and the resolvent of the adjacency tensor
associated with this network using the Einstein product. We showed how to
approximate these tensor functions via some Krylov subspace methods based on
tensor format. Numerical tests lead to satisfactory results. Exploring general
multilayer networks using tensorial formalism to analyze and extract more
information about this kind of networks is a future project.

## References

  * [1] M. Al Mugahwi, O. De la Cruz Cabrera, and L. Reichel, Orthogonal expansion of network functions, Vietnam J. Math., 48 (2020), pp. 941–962.
  * [2] B. Beckermann and L. Reichel, Error estimation and evaluation of matrix functions via the Faber transform, SIAM J. Numer. Anal., 47 (2009), pp. 3849–3883.
  * [3] F. P. A. Beik, K. Jbilou, M. Najafi-Kalyani, and L. Reichel, Golub-Kahan bidiagonalization for ill-conditioned tensor equations with applications, Numer. Algorithms, 84 (2020), pp. 1535–1563.
  * [4] A. Bentbib, M. El Ghomari, K. Jbilou, and L. Reichel, The Golub-Kahan method and Gauss quadrature for tensor function approximation, Numer. Algorithms, in press.
  * [5] M. Benzi and C. Klymko, Total communicability as a centrality measure, J. Complex Netw., 1 (2013), pp. 124–149.
  * [6] K. Bergermann, Multiplex-matrix-function-centralities,
https://github.com/KBergermann/Multiplex-matrix-function-centralities.

  * [7] K. Bergermann and M. Stoll, Fast computation of matrix function-based centrality measures for layer-coupled multiplex networks, Phys. Rev. E, 105 (2022), Art. 034305.
  * [8] R. Behera and D. Mishra, Further results on generalized inverses of tensors via the Einstein product, Linear Multilinear Algebra, 65 (2017), pp. 1662–1682.
  * [9] M. Brazell, N. Li, C. Navasca, and C. Tamon, Solving multilinear systems via tensor inversion, SIAM J. Matrix Anal. Appl., 34 (2013), pp. 542–570.
  * [10] S. Cipolla, M. Redivo-Zaglia, and F. Tudisco, Shifted and extrapolated power methods for tensor $\ell^{p}$-eigenpairs, Electron. Trans. Numer. Anal., 53 (2020), pp. 1–27.
  * [11] M. De Domenico, A. Solé-Ribalta, E. Cozzo, M. Kivelä, Y. Moreno, M. A. Porter, S. Gómez, and A. Arenas, Mathematical formulation of multilayer networks, Phys. Rev. X, 3 (2013), Art. 041022.
  * [12] M. De Domenico, A. Solé-Ribalta, E. Omodei, S. Gómez, and A. Arenas, Centrality in interconnected multilayer networks, arXiv:1311.2906v1, (2013).
  * [13] O. De la Cruz Cabrera, M. Matar, and L. Reichel, Analysis of directed networks via the matrix exponential, J. Comput. Appl. Math., 355 (2019), pp. 182–192.
  * [14] M. El Guide, A. El Ichi, F. P. Beik, and K. Jbilou, Tensor Krylov subspace methods via the Einstein product with applications to image and video processing, Appl. Numer. Math., 181 (2022), pp. 347–363.
  * [15] S. El-Haloy, S. Noschese, and L. Reichel, Perron communicability and sensitivity of multilayer networks, Numer. Algorithms, in press.
https://doi.org/10.1007/s11075-022-01418-4

  * [16] E. Estrada, Communicability geometry of multiplexes, New J. Phys., 21 (2019), Art. 015004.
  * [17] E. Estrada, The Structure of Complex Networks: Theory and Applications, Oxford University Press, Oxford, 2011.
  * [18] E. Estrada and N. Hatano, Communicability in complex networks, Phys. Rev. E, 77 (2008), Art. 036111.
  * [19] E. Estrada and D. J. Higham, Network properties revealed through matrix functions, SIAM Rev., 52 (2010), pp. 696–714.
  * [20] E. Estrada and J. A. Rodriguez-Velazquez, Subgraph centrality in complex networks, Phys. Rev. E, 71 (2005), Art. 056103.
  * [21] C. Fenu, D. Martin, L. Reichel, and G. Rodriguez, Block Gauss and anti-Gauss quadrature with application to networks, SIAM J. Matrix Anal. Appl., 34 (2013), pp. 1655–1684.
  * [22] G. H. Golub and G. Meurant, Matrices, Moments and Quadrature with Applications, Princeton University Press, Princeton, 2010.
  * [23] B. Huang, X. Yajun, and C. Ma, Krylov subspace methods to solve a class of tensor equations via the Einstein product, Numer. Linear Algebra Appl., 26 (2019), Art. e2254.
  * [24] K. Jbilou, A. Messaoudi, and H. Sadok, Global FOM and GMRES algorithms for matrix equations, Appl. Numer. Math., 31 (1999), pp. 49–63.
  * [25] K. Jbilou, H. Sadok, and A. Tinzefte, Oblique projection methods for multiple linear systems, Electron. Trans. Numer. Anal., 20 (2005), pp. 119–138.
  * [26] L. Katz. A new status index derived from sociometric analysis, Psychometrika, 18 (1953), pp. 39–43.
  * [27] T. G. Kolda and B. W. Bader, MATLAB tensor toolbox, Sandia National Laboratories (SNL), Albuquerque, NM, and Livermore, CA, 2006.
  * [28] T. G. Kolda and B. W. Bader, Tensor decompositions and applications, SIAM Rev., 51 (2009), pp. 455–500.
  * [29] M. Kivelä, A. Arenas, M. Barthelemy, J. P. Gleeson, Y. Moreno, and M. A. Porter, Multilayer networks, J. Complex Netw., 2 (2014), pp. 203–271.
  * [30] K. Lund, The tensor t-function: A definition for functions of third-order tensors, Numer. Linear Algebra Appl., 27 (2020), Art. e2288.
  * [31] L. Qi and Z. Luo, Tensor Analysis: Spectral Theory and Special Tensors, SIAM, Philadelphia, 2017.
  * [32] D. Taylor, M. A. Porter, and P. J. Mucha, Tunable eigenvector-based centralities for multiplex and temporal networks, Multiscale Model. Simul., 19 (2021), pp. 113–147.
  * [33] M. Wu, S. He, Y. Zhang, J. Chen, Y. Sun, Y. Liu, J. Zhang, and H. V. Poor, A tensor-based framework for studying eigenvector multicentrality in multilayer networks, Proc. NAS, 116 (2019), pp. 15407–15413.


It would be ideal to create a sound and complete procedure to verify that any quantum program is synchronized. In the theoretical worst case, this problem has exponential complexity, as is the case for verifying typical nontrivial semantic properties of classical programs with bounded registers.
Nevertheless, for practical instances of programs for the quantum control machine, it is unnecessary to determine the results of an exponential number of executions.

First, verifying a sequence of $n$ binary conditional branches requires checking only that each of the $n$ pairs of syntactic branches align, rather than the $2^n$ possible program paths.
Second, we present in \Cref{sec:iteration-case-study} a reasoning strategy for loops in which the execution time depends on only a classical upper bound, making the program synchronized by construction. Constant-time compilation strategies in secure programming~\citep{barthe2020} similarly apply.

In a more general setting, we expect that given sufficient analysis of program structure, it will be possible to partially decide this condition through techniques developed in classical and quantum program analysis. Five promising approaches are:

\begin{enumerate}
\item Formally verify, only once, the compilation of high-level constructs such as branches and loops, similarly to verified optimizers for quantum circuits~\citep{hietala2021};
\item Generalize the technique of abstracting over quantum variables as in the loop case above, similarly to quantum abstract interpretation~\citep{yu2021};
\item Use the interpretation of synchronization as a $k$-safety property to leverage relational and Cartesian Hoare logics that provide automated verification tools~\citep{sousa2016}, as could be done via quantum relational Hoare logic~\citep{unruh2019};
\item Use the interpretation of synchronization as a bound on the worst-case execution time of a program to leverage analyses of execution time via SMT solvers~\citep{seshia2011};
\item Perform runtime checks during execution on quantum hardware using a small number of measurements~\citep{li2020,twist}.
\end{enumerate}

\section{Other Related Work} \label{sec:related-work}

In this section, we discuss related work that was not addressed in \Cref{sec:implications}, on the topics of quantum $\lambda$-calculi, reversible and oblivious computation, and unbounded-time quantum computation.

\paragraph{Quantum $\lambda$-Calculi}
Control flow in superposition hints at the possibility of a purely functional paradigm of quantum programming.
Though several $\lambda$-calculi for quantum computation have been proposed, their operational semantics are not inherently injective, and they do not permit $\lambda$-terms to exist in superposition.
For example, the quantum $\lambda$-calculus with classical control~\citep{selinger2004} performs $\beta$-reduction classically and does not permit terms to exist in superposition.
Another example is the calculus of \citet{sabry2018}, whose operational semantics of $\lambda$-terms is not injective as both $\lambda x.x$ and $(\lambda x.x)(\lambda x.x)$ reduce to $\lambda x.x$. A possible way toward developing an inherently injective calculus may be to use a type system to reject terms that would reduce in a non-injective manner, as has been pursued by~\citet{arrighi2017,diazcaro2019}.

To our knowledge, the understanding of higher-order quantum programming by categorical approaches~\citep{pagani2014,malherbe2013,rennela2017,hines2011} remains that ``the general connection between physics and higher-order quantum functions is yet unclear''~\citep{renella2018}. Recent developments~\citep{hasuo2017,clairambault19,clairambault20,kornell2021} may be able to shed further light. A related topic is quantum computation with indefinite causal order such as the quantum switch~\citep{chiribella2013}, for which the quantum control machine could serve as a programmable instruction set.

\paragraph{Reversible Computation}
In classical reversible computation, the Landauer embedding~\citep{landauer1961} and uncomputation~\citep{bennett1973} have long been used to simulate arbitrary control flow constructs, including rewriting systems such as combinatory logic~\citep{dipierro2006}.

By contrast, the no-embedding theorem implies that they fail to correctly do so in quantum computation. The reason is entanglement, an inherently quantum phenomenon that does not arise in reversible computation.
The theorem precludes any system of uncomputation that is program-agnostic and always correctly erases the history in the Landauer embedding. It does not preclude program-specialized or sometimes-incorrect systems, and indeed reverse jumps are instances of manual uncomputation~\citep{silq,tower} in quantum programming.

The branch control register was introduced by seminal work on the Pendulum reversible architecture~\citep{axelsen2007,vieri1998,frank1999,axelsen2011}. In reversible computation, the cost of using a history is extra memory or energy consumption~\citep{thomsen2012}. By contrast, in quantum computation, entanglement leads the output to be irrecoverably incorrect.

\paragraph{Oblivious Computation}
Synchronization is related to the property of \emph{obliviousness}~\citep{pippenger1979,goldreich1996} as studied in security and privacy. Obliviousness states that the evolution of the program counter may not depend on the input data, only on its size, and is used in the proof~\citep[Theorem 4.3]{bv1997} that any quantum Turing machine may be transformed into an equivalent one with a classically fixed termination time.

Compared to synchronization, obliviousness is unnecessarily strong --- it constrains all intermediate states of the computation, meaning that the program in \Cref{sec:examples} could not manipulate the program counter based on data if it were required to be oblivious. Synchronization thus lies between the properties of obliviousness and \emph{history independence}~\citep{teague2001,tower}, which fixes the final states, but not running times, of a computation.

\paragraph{Unbounded-Time Quantum Computation}
A restriction of the quantum control machine is that the running time of a computation must be a classically determined value. This restriction is fundamentally implied by synchronization, which states that there must exist a point in time at which all branches of the machine state superposition possess the same control state.

Evidence suggests that such a restriction is inevitable to quantum computation, and that it is not possible to physically realize a quantum loop that does not depend on any classical bound:

\begin{itemize}
\item \citet{myers1997,kieu1998,linden1998} show that in a quantum Turing machine, when different branches of the state superposition of a computation take different amounts of time to terminate, no general scheme is known for determining whether all branches have halted that does not collapse the superposition.
\item \citet{ying2010} indicate that ``If a quantum loop is allowed to be embedded into another quantum loop, then \ldots{} the body of the latter loop is not a unitary operator but a super-operator in general'', making it not realizable as a circuit of unitary logic gates.
\item \citet{andresmartinez2022} show that a loop-like operation is possible via weak measurements of the system. The issue is that if the strength parameter used by the scheme is too strong or weak, the result is incorrect, and in turn determining its optimal value requires access to problem structure analogous to a classical iteration bound.
\end{itemize}

As a result, the quantum control machine has been designed to accord with the standard assumption in the theory of quantum computation~\citep[Section 3.5]{bv1997} that the termination time of a quantum computation is classically fixed.

\end{document}
\endinput