Faster Tensor Canonicalization

TL;DR

This paper introduces a modified tensor canonicalization algorithm that significantly improves performance for common symmetric or antisymmetric tensor expressions, reducing complexity from factorial to polynomial time in typical cases.

Contribution

A new algorithm and implementation that optimize tensor canonicalization for symmetric cases, reducing computational complexity and increasing speed in practical scenarios.

Findings

Polynomial time performance for common symmetric tensor cases

Significant speed improvements over previous algorithms

Limited worst-case factorial complexity in rare scenarios

Abstract

The Butler-Portugal algorithm for obtaining the canonical form of a tensor expression with respect to slot symmetries and dummy-index renaming suffers, in certain cases with a high degree of symmetry, from $O(n!)$ explosion in both computation time and memory. We present a modified algorithm which alleviates this problem in the most common cases---tensor expressions with subsets of indices which are totally symmetric or totally antisymmetric---in polynomial time. We also present an implementation of the label-renaming mechanism which improves upon that of the original Butler-Portugal algorithm, thus providing a significant speed increase for the average case as well as the highly-symmetric special case. The worst-case behavior remains $O(n!)$, although it occurs in more limited situations unlikely to appear in actual computations. We comment on possible strategies to take if the nature of a computation should make these situations more likely.

Tables

Table 1. Table 1 : Group codes for Label-Groups array.

Code	Meaning
$\varnothing$	No label exchange
$\sim$	Component labels such as $⟨1_1,1_2,1_3,1_4⟩$
$+$	Dummy labels with symmetric metric such as $g_{ab}$
$-$	Dummy labels with antisymmetric metric such as ${\epsilon }_{\alpha \beta }$
$\perp$	Lower dummy labels with no metric
$\top$	Upper dummy labels with no metric

Equations

$$\mathbf{v}=v^i{\mathbf{e}}_i\equiv v^1{\mathbf{e}}_1+v^2{\mathbf{e}}_2+\dots +v^n{\mathbf{e}}_n,$$

$$T^{\circ }{}_{\circ \circ \circ },$$

$$(M^{\top }{)}_{ij}\equiv M_{ji}.$$

$$T^1{}_{223}$$

$$T_{ab}{}^{ab}=T_{bc}{}^{bc}.$$

$$R_{abcd}=-R_{bacd},R_{abcd}=R_{cdab},$$

$$R_{abcd}+R_{bcad}+R_{cabd}=0.$$

$$T_{ab}{}^{ab}=T_{ba}{}^{ba},$$

$$A_{11}=A_{22}=0,A_{12}=-A_{21}.$$

$$T_{a_1{a}_2\dots a_n}=\epsilon T_{\pi (a_1{a}_2\dots a_n)}.$$

$$T_{a_1{a}_2\dots a_n}-{\epsilon }_1{T}_{{\pi }_1(a_1{a}_2\dots a_n)}-{\epsilon }_2{T}_{{\pi }_2(a_1{a}_2\dots a_n)}-\dots -{\epsilon }_k{T}_{{\pi }_k(a_1{a}_2\dots a_n)}=0.$$

$$T_{bdafce}$$

$$g=\left(\begin{array}{cccccc}1& 2& 3& 4& 5& 6\\ b& d& a& f& c& e\end{array}\right).$$

$$g^{\prime }=\left(\begin{array}{cccccc}3& 1& 2& 4& 5& 6\\ b& d& a& f& c& e\end{array}\right)=\left(\begin{array}{cccccc}1& 2& 3& 4& 5& 6\\ d& a& b& f& c& e\end{array}\right)\Rightarrow T_{dabfce}.$$

$$T_{11ab}{}^{bc}=T_{1_1{1}_{2}a{b}_1}{}^{b_{2}c},$$

$$\left(\begin{array}{cccccc}1& 2& 3& 4& 5& 6\\ 1_1& 1_2& a& b_1& b_2& c\end{array}\right).$$

$$\left(\begin{array}{cccccc}1& 2& 3& 4& 5& 6\\ 1_2& 1_1& a& b_1& b_2& c\end{array}\right),\left(\begin{array}{cccccc}1& 2& 3& 4& 5& 6\\ 1_1& 1_2& a& b_2& b_1& c\end{array}\right),\left(\begin{array}{cccccc}1& 2& 3& 4& 5& 6\\ 1_2& 1_1& a& b_2& b_1& c\end{array}\right).$$

$${\psi }^{\alpha }{\chi }_{\alpha }=-{\psi }_{\alpha }{\chi }^{\alpha }.$$

$$g=\left(\begin{array}{cccccccc}1& 2& 3& 4& 5& 6& 7& 8\\ b& d& a& f& c& e& +& -\end{array}\right),$$

$$g^{\prime }=\left(\begin{array}{cccccccc}1& 2& 3& 4& 5& 6& 7& 8\\ d& b& a& f& c& e& -& +\end{array}\right),$$

$$g=⟨b,d,a,f,c,e,+,-⟩,$$

$$g^{\prime }=\ell [[g[[s]]]],$$

$$g[[s]][[i]]=g[[s[[i]]]]=i^{s\cdot g}.$$

$$T_a{}^a{}_b{}^b{}_c{}^c,$$

$$L=\{+(1,2),+(3,4),+(5,6),+(1,3)(2,4),+(3,5)(4,6)\},$$

$$L=\{+(3,4),+(1,2),+(5,6),+(1,3)(2,4),+(1,5)(2,6)\}.$$

$$T_{1111},$$

$$L=\{+(1,2),+(2,3),+(3,4)\},$$

$$T_{abc}=T_{bca}=T_{cab},$$

$$T_{cb}{}^c{U}^{ba}{}_a,$$

Full text

1. 1 Introduction

2. 2 Tensor symmetries and permutation groups

3. 2.1 Specification of the problem

4. 3 The Butler-Portugal algorithm

5. 3.1 The label symmetry group

6. 3.2 The slot symmetry group

7. 3.3 The algorithm

8. 3.4 Worst-case complexity analysis

9. 4 The improved algorithm

10. 4.1 Insights from Penrose graphical notation

11. 4.2 Data structures

12. 4.3 The main algorithm

13. 4.4 Comments on complexity

14. 5 Performance testing

15. 5.1 Preliminaries

16. 5.2 Random Riemann invariants

17. 5.3 Total symmetry

18. 5.4 Pairwise symmetry

19. 6 Discussion

20. A Efficient computation with permutation groups

21. B Subprocedures of the improved algorithm

1 Introduction

Computer algebra systems have become quite powerful at simplifying expressions with polynomials, integrals, derivatives, special functions, etc. All of these fall under the category of what we will call scalar algebra, meaning that they deal with objects that are functions from $\mathbb{R}\to\mathbb{R}$ or $\mathbb{C}\to\mathbb{C}$. A field that is still quite in development is the extension of this computational power to tensor algebra, by which we mean the manipulation of objects with “indices” such as vectors $v^{i}$, matrices $M_{ij}$, and higher-rank tensors $T^{i}{}_{jk\ell}$. Such indexed objects are a way to denote multidimensional arrays (and thus represent objects of multilinear algebra). The indices refer to the components of the object in a basis of some vector space like $\mathbb{R}^{n}$; for example, one can write:

[TABLE]

where $\mathbf{v}\in\mathbb{R}^{n}$ is a vector and $\{{\mathbf{e}_{i}}\}$ is a given basis.111 “Upstairs” (superscript) and “downstairs” (subscript) indices are to be distinguished, and when the same index is repeated both upstairs and downstairs, it is meant to be contracted or summed over, as shown in (1.1). Where multiple indices occur on a single tensor, these refer to the component along the tensor product of the corresponding basis elements, as in $\mathbf{T}=T^{i}{}_{jk\ell}\,\mathbf{e}_{i}\otimes\mathbf{e}^{j}\otimes\mathbf{e}^{k}\otimes\mathbf{e}^{\ell}$. Alternatively, we can choose to think of indices abstractly, referring merely to the structure of a tensor and its contractions, rather than literally referring to the components of the underlying multidimensional array [1].In either case, whether indices are thought of as literal or abstract, the algebra of indexed objects is as follows: We can think of each index as a slot which has been filled with a label such as $i,j,k,\ell$, etc. The slots are, for our purposes, an ordered sequence of “blanks” in a tensor expression into which labels can go,

[TABLE]

whereas the labels $i,j,k,\ell$ are a means to “name” these slots, and either indicate that a pair of slots have been contracted (in the case that a label is repeated), or indicate that a slot is free; that is, available to be contracted.222Thus in $T^{k}{}_{ikj}$, the first and third slots are contracted, while the second and fourth slots are free. We will frequently refer to the slots by their names, saying that $i,j$ are free, while $k$ is contracted. Note that a given label may only be repeated twice, once upstairs and once downstairs, to indicate a single contraction. The positional order of the labels has meaning, and can be used to effectively indicate notions like the transpose of matrices:

[TABLE]

The free labels must be “balanced” in any sum or equation in order for it to be a valid mathematical sentence. For repeated labels which indicate a contraction, one may use any available symbol; one calls these “dummy” labels, as the particular symbol used is unimportant. A final use of labels is one we will call component labels, which refer to particular components of a tensor (in a particular concrete basis); for example,

[TABLE]

refers to the component of $\mathbf{T}$ along the tensor product basis element $\mathbf{e}_{1}\otimes\mathbf{e}^{2}\otimes\mathbf{e}^{2}\otimes\mathbf{e}^{3}$.The task of a tensor computer algebra system—that is, to manipulate algebraic expressions with indexed objects—is complicated by two types of symmetries: re-labelling symmetries, which we have already hinted at, and intrinsic symmetries, which we will define in a moment. Re-labelling symmetries arise when there are dummy labels present, as illustrated by

[TABLE]

Since the labels used to indicate each contraction are unimportant, the two sides of (1.5) are equal, even though they do not look the same. In order to see that they are equal, the system must either canonicalize the sequence of dummy labels used (thus putting each monomial into a standard form), or use some sort of isomorphism-detection algorithm (in order to match one monomial to the other).With only dummy-relabelling to contend with, this problem would be rather simple. However, a tensor can also have intrinsic symmetries, which are symmetries that involve permuting the tensor’s slots. For example, a symmetric matrix $M_{ij}=M_{ji}$ is symmetric irrespective of the particular labels used; one can think of this symmetry as a property of $\mathbf{M}$ itself. Intrinsic symmetries come in two distinct forms, which we will call mono-term symmetries and multi-term symmetries. Some simple examples of mono-term symmetries are the slot-exchange symmetries of the Riemann curvature tensor,

[TABLE]

whereas an example of a multi-term symmetry is the algebraic Bianchi identity:

[TABLE]

In the context of a calculation (that is, when labels have been placed in the slots), the intrinsic symmetries of a tensor monomial can interact with its re-labelling symmetries in a complicated way. For example, a tensor $T_{abcd}$ which has no intrinsic symmetries may acquire the appearance of a slot symmetry when contracted,

[TABLE]

because we are free to exchange dummy labels. Similarly, if a tensor $A_{ab}=-A_{ba}$ is intrinsically antisymmetric, then placing component labels into its slots can have non-trivial results:

[TABLE]

Here we can think of the first equalities as arising from the contradiction between the antisymmetry of $A$ and the symmetry of exchanging identical component labels. The same sort of situation can happen with dummy labels, if say we take $M^{ab}A_{ab}$ where $M$ is symmetric and $A$ is antisymmetric. In general, all of these complications may be mixed within a single expression: a tensor monomial may have intrinsic slot symmetries (of both single- and multi-term types), some free indices, some dummy indices, and some concrete component labels, all at the same time.Thus any CAS with the power to simplify tensor expressions is faced with a challenge to efficiently manipulate index symmetries. The greatest amount of progress has been made on mono-term symmetries, whereas multi-term symmetries are considerably more difficult to manage. For mono-term symmetries, there are two main approaches: index canonicalization and index isomorphism:The index canonicalization approach seeks to reduce each monomial in an expression to its canonical form, which is defined as the least possible arrangement of its indices, modulo slot symmetries and label exchanges, where least is defined by lexicographical order. Simplification can then proceed as it does in a standard (i.e. non-tensorial) CAS, as all objects which are equal have been reduced to the same form. While processing a given monomial, the same algorithm can easily detect if the monomial is itself equal to zero due to its own internal symmetries (say, by having a pair of symmetric slots which are contracted against a pair of antisymmetric ones); thus, the canonicalization approach can also pre-emptively simplify some parts of expressions. The canonicalization strategy is well-suited for use as a plugin or package for a standard CAS such as Maple or Mathematica, since it reduces terms to a form which can then be manipulated by pre-existing routines.By contrast, the index isomorphism approach works by attempting to match the various tensor monomials of an expression onto each other by constructing a bijection between their free indices (subject to symmetries). With every successful bijection, the expression can be simplified by combining the corresponding terms. This approach has two advantages: First, the algorithm to make a single match could conceivably run faster than a canonicalization algorithm, since one has a definite goal in mind rather than needing to search the space of equivalent configurations for the “least” one.333At the time of writing, however, we do not know of any software which actually implements such an algorithm which is efficient in the presence of slot symmetries. Second, a matching algorithm can treat compound expressions such as $T^{ac}(S^{b}{}_{bcd}+V_{cd})$ as a single monomial $(\cdot)^{a}{}_{d}$, and thus simplify the expression tree without flattening it; whereas the canonical form of such an expression is not obvious, and may be ill-defined.444However, the notion of canonical-ness of compound expressions has not been explored in the literature, nor implemented in any software package to our knowledge, and it may well be that one can come up with a workable definition! However, the disadvantages of the isomorphism approach may outweigh the advantages: First, while a single matching test might run faster than canonicalization, one still has to make up to $O(n^{2})$ matching tests, vs. $O(n)$ canonicalizations (where $n$ here is the number of terms in an expression). Second, a matching algorithm alone has no good way of detecting whether a single monomial vanishes due to its own symmetry, so in the end one has to run a separate test for that. And finally, attempting to use the index isomorphism approach within an existing (scalar) CAS is cumbersome and requires circumventing what that system already knows about expression simplification; it is an approach better suited to standalone software.There are many implementations of the index canonicalization strategy in currently-available CAS software.555There is currently only one software package which uses the index isomorphism strategy: the standalone tensor CAS Redberry [8]. Some of the more popular ones include the open-source Canon package [2] for Maple, which is also included in the GRtensorII [3] suite of packages; the open-source xPerm package [4] of the package suite xAct [5] for Mathematica; the standalone open-source tensor CAS Cadabra [6] (which in fact uses the same code from xPerm); and the tensor subsystem of the standalone open-source Python CAS SymPy [7]. Since Version 9, Mathematica also includes a built-in implementation of tensor canonicalization, which is not open-source, but we can nevertheless draw conclusions about its algorithm from testing.The most efficient algorithm in the literature for index canonicalization is the Butler-Portugal algorithm [9, 10], which is used in all of the above-mentioned software packages, and has been shown to have polynomial complexity (in the number of indices) in most common cases. However, in tensor calculus in higher dimensions (of particular interest to the high-energy theoretical physics community), there are certain types of simplification problems which thwart the Butler-Portugal algorithm and cause it to consume $O(n!)$ memory and take $O(n!)$ time. These are problems where tensors with many indices and a high degree of symmetry are contracted with each other; thus the slots of these highly-symmetric tensors are filled with dummy indices and are subject to a high degree of re-labelling ambiguity. The overlap of slot symmetries with re-labelling symmetries inhibits the algorithm’s ability to deduce whether two arrangements of indices are actually different.In this article we present a modification which improves upon the Butler-Portugal algorithm by automatically recognizing (in polynomial time) these worst-case situations, and adapting to avoid processing $O(n!)$ redundant steps. Thus for tensor monomials with total symmetry or total antisymmetry over some subset(s) of slots, our improved algorithm can canonicalize them in polynomial time, even subject to index re-labelling symmetries.666We note that some attempt has been made at improving the canonicalization of (anti)symmetric tensors in [11], although it is far from complete. In the course of implementing these modifications for these worst-case situations, we introduce a number of data structures which significantly increase the algorithm’s average-case performance as well.Aside from (subsets of indices with) total symmetry or total antisymmetry, there are still situations that can lead to $O(n!)$ behavior (such as a symmetric group which acts on pairs of indices rather than individuals), and there is no efficient means to identify all potential sources of combinatorial explosion. Thus our algorithm still has worst-case time and space complexities of $O(n!)$, although the cases in which this behavior occurs are much less frequent in practical calculations. However, in the event that one does expect these situations to occur, we discuss some possible strategies based upon the concepts introduced in this work.This paper is structured as follows: in Section 2 we review tensor symmetries and permutation groups, and precisely define the problem to be solved; in Section 3, we review the Butler-Portugal algorithm; in Section 4 we present the improved algorithm; in Section 5 we give the results of performance testing on a few different types of tensor inputs; and in Section 6 we discuss these results and potential further improvements.

2 Tensor symmetries and permutation groups

In this work we will focus on the mono-term symmetries of indexed objects and their relation to permutation groups. A mono-term symmetry of a tensor can be described by a signed permutation777One could also consider promoting $\varepsilon$ to be a root of unity. $(\varepsilon,\pi)\in\{{-1,1}\}\times S_{n}$, where $S_{n}$ is the symmetric group on the tensor’s slots:

[TABLE]

That is, the tensor returns to itself, up to a sign, after permuting its slots in some specific way. This type of symmetry encompasses symmetric and antisymmetric matrices, and most of the examples discussed in the Introduction.888By contrast, to describe a multi-term symmetry of a tensor (such as the one in the algebraic Bianchi identity (1.7)) requires several signed permutations $\{{(\varepsilon_{1},\pi_{1}),(\varepsilon_{2},\pi_{2}),\ldots,(\varepsilon_{k},\pi_{k})}\}$, and takes the form

$T_{a_{1}a_{2}\ldots a_{n}}-\varepsilon_{1}\,T_{\pi_{1}(a_{1}a_{2}\ldots a_{n})}-\varepsilon_{2}\,T_{\pi_{2}(a_{1}a_{2}\ldots a_{n})}-\ldots-\varepsilon_{k}\,T_{\pi_{k}(a_{1}a_{2}\ldots a_{n})}=0.$

Written in this way, we see that (2.1) arises as a special case. Dealing with multi-term symmetries efficiently is an important open problem in the symbolic computation literature, which we will not address here. One approach used in the Cadabra system [6] is to enforce the multi-term symmetry by projection onto Young tableaux, but this has the tendency to introduce many extra terms into an expression. Another approach is described in [12] and is based on graph isomorphism.

2.1 Specification of the problem

For mono-term symmetries, one can address the canonicalization problem purely from the standpoint of permutation groups [10], which are well-studied [9, 13, 14]. Thus using this well-honed tool, we will focus on the mono-term symmetries.First let us formalize the definition of the problem. It is useful first of all to think of a tensor monomial

[TABLE]

as a map from a space of slots into a space of labels, which therefore informs us which label should go into which slot. We will refer to this map as a configuration $g:\mathcal{S}\to\mathcal{L}$, where $\mathcal{S}$ is the space of slots and $\mathcal{L}$ is the space of labels. The configuration corresponding to (2.2) can then be written as

[TABLE]

We note that $g$ is essentially a permutation—if one instead replaces the labels $a$ through $f$ by numbers 1 to 6, then the notation (2.3) becomes the standard Cauchy two-line notation for a permutation map [15]. However, we will find it useful to maintain a distinction between the slot and label spaces.If a tensor monomial has mono-term symmetries which permute the slots, we can think of these as acting on the first row of $g$. Thus if we act first with a permutation (in cyclic notation) to exchange slots $(1,2)$, and then afterwards act with one to exchange $(2,3)$, we obtain

[TABLE]

The slot symmetries of a tensor monomial form a group $S$ which acts on the slot space $\mathcal{S}$.If a tensor monomial has multiple repeated labels, such as when there are dummy pairs or concrete component labels, we may clarify what happens by treating the repeated labels as distinct. We can then implement their “sameness” via symmetries which exchange the labels. For example, we could have

[TABLE]

where we attach a subscript to some of the labels to distinguish identical ones999For dummy pairs we adopt the convention that the “downstairs” label comes first. This is opposite the convention of [10, 4].. Then the configuration could be given as

[TABLE]

Due to the repetition of labels, even if $T$ has no slot symmetries, the following configurations are equivalent to (2.6) (assuming the presence of a metric tensor):

[TABLE]

Thus we see that repeated labels of either type give us a symmetry group $L$ that acts on label space $\mathcal{L}$. These symmetries are not intrinsic to the tensor, but arise from the context; i.e., having placed labels in the slots of the tensor.The mono-term symmetries of a tensor monomial then consist of the union of two distinct symmetry actions: those symmetries, intrinsic to the definition of $T$, which act on the slots $\mathcal{S}$; and those symmetries, arising from context, which act on the labels $\mathcal{L}$. Again treating the configuration $g$ as a permutation, the group of slot symmetries $S$ act from the right of $g$, and the group of label symmetries $L$ act from the left. Thus the set of all equivalent configurations of the indices of a tensor monomial is the double coset $LgS$, and the task of index canonicalization is precisely that of finding the canonical representative of a double coset.In all of this discussion thus far, we have ignored the fact that tensor symmetries are signed permutations. In fact, both groups $S$ and $L$ can have signed permutations; the ones in $L$ can come from “antisymmetric metrics”, used in spinor formalisms in certain dimensions, which can introduce minus signs when raising/lowering contracted indices, e.g.

[TABLE]

This corresponds, in label space, to the signed permutation $-(\alpha_{1},\alpha_{2})$ that exchanges the two (distinguished) $\alpha$’s. It is not difficult to expand our formalism to include signed permutations. As was observed in [4], the sign is just an element of a $\mathbb{Z}_{2}$ group, which can be represented as a permutation group. One just needs to add an extra two columns at the end of $g$:

[TABLE]

and in label space, allocate two arbitrary symbols for tracking the order of the last two columns. The action of the signed permutation $-(1,2)$ on (2.9) then gives

[TABLE]

and thus, one can work out whether the sign of any expression has flipped by checking the last pair of columns.We point out that $g$ has a fairly obvious computer implementation which makes clear our notions of left- and right-actions. If we store the configuration of (2.3), (2.9) in computer memory as the array101010Of course, to encode the symbols of the label alphabet, we should use numbers in the same range as the slot numbers, so this becomes $g=\langle{2,4,1,6,3,5,7,8}\rangle$. So in the end, $g$ is truly just a permutation map.

[TABLE]

then the map $g$ can be implemented as array access $i\mapsto g\llbracket i\rrbracket$, for $i\in\mathcal{S}$. One can easily extend this notion to the group multiplication of permutations, and thus for $\ell\in L$ and $s\in S$, one has schematically

[TABLE]

which corresponds to the stated notions of left and right.111111Alternatively, if one wants to use cyclic notation for permutations, it is more convenient to work with the right-action $i\mapsto i^{g}$, in which case the notions of left and right are reversed: $g^{\prime}=s\cdot g\cdot\ell$. Finally, one of the main advantages of treating the configuration as a slots-to-labels map, rather than labels-to-slots, is that it becomes trivial to order configurations lexicographically; for example one has that $\langle{b,a,f,c,e,d,+,-}\rangle\prec\langle{b,c,a,d,f,e,+,-}\rangle$, etc., for the lexicographical ordering $\prec$. Our choice to append the sign columns to the end rather than the beginning is motived by this lexicographic ordering: if a set of configurations contains both $+g$ and $-g$, then these two configurations will be adjacent when the set is sorted lexicographically. Thus it becomes easy to detect when the only tensor consistent with a given set of index permutations is the zero tensor.

3 The Butler-Portugal algorithm

In [10], Manssur and Portugal adapt Butler’s algorithm [9] for finding the canonical representative of a double coset to the specific problem of simplifying tensor monomials with dummy indices. Their main modification, which results in a highly-efficient algorithm for most tensor monomials in practice, is to take advantage of the simple and predictable nature of the label permutation group $L$ for dummy indices. In [4], Martín-García extends the algorithm of [10] to include dummy indices drawn from multiple separate pools (in the case that a tensor has indices belonging to different vector bundles), as well as component labels. The label permutation group $L$ retains a simple structure, which allows these modifications to run equally efficiently.Like many algorithms in finite group theory, the Butler-Portugal algorithm makes use of the concepts of a base and strong generating set to deal with the permutation groups in an efficient way. We refer the reader to Appendix A for a brief presentation of these concepts.

Notational conventions

Before delving into algorithms, we briefly mention some notational conventions. We will use curly braces $\{{\ldots}\}$ to denote sets, in the strict sense of being unordered. To denote arrays or lists, whose elements are ordered by their position, we use angle brackets $\langle{\ldots}\rangle$. To denote array access, we use double square brackets; thus $g\llbracket i\rrbracket$ is the $i$-th entry of $g$. Double square brackets can also denote the action of a permutation, if their argument is itself a permutation map rather than a single number; thus $g\llbracket s\rrbracket$ denotes the permutation that results by acting with $g$ on $s$. The convention for $\llbracket\cdot\rrbracket$ is left multiplication; whereas in cyclic notation, the equivalent action is right multiplication. Thus the array $g\llbracket s\rrbracket$ is the same as the array corresponding to $s\cdot g$ in cyclic notation; that is the action on some element $i$ is “first $s$, then $g$”:

[TABLE]

In the pseudocode of the algorithms we present, we will adopt a few conventions on typefaces. Variables will be indicated by sans-serif type, and arrays will be indicated by Sans-Serif-with-Caps. Function names will take Small-Caps-Camel-Case.Finally, we will frequently refer to “lexicographical ordering”, by which we mean an ordering $\prec$ in which all of the free indices come first, and then the component indices and dummy indices of various types. For dummy indices, the ordering is defined such that they come in pairs, lower first and then upper. Therefore the indices of the tensor $T^{ab}{}_{fgc}{}^{c}{}_{d}{}^{d}{}_{e}{}^{e}$ are in lexicographical order.

3.1 The label symmetry group

The observation of [10] which contributes most to the speed of its algorithm is that the label symmetry group $L$ takes a simple form. As an example, take the tensor monomial

[TABLE]

whose indices come only in contracted pairs. Acting on the label space $\mathcal{L}=\{{a_{1},a_{2},b_{1},b_{2},c_{1},c_{2}}\}$, we can immediately see that the label symmetry group is given by the signed permutations

[TABLE]

and moreover, this is a strong generating set with respect to the base $B=\langle{1,3,5}\rangle$. More importantly, however, it is simple to re-order the base: If instead we want a strong generating set with respect to the base $B=\langle{3,1,5}\rangle$, then we need only conjugate the set (3.3) by the signed permutation $+(1,3)(2,4)$, which obtains

[TABLE]

Butler’s algorithm [9] for the canonical representative of a double coset relies on repeatedly changing the base order for one of the two subgroups. The best known algorithms for an arbitrary base rearrangement on generic groups require $O(n^{3})$ time [16, 17]. The conjugation operation, however, is much faster. Portugal observed [10] that for tensor canonicalization, the fixed, simple structure of the label symmetry group always allows for a much-streamlined base-change-by-conjugation.A similar situation occurs with component labels. This time take the tensor monomial

[TABLE]

with label space $\mathcal{L}=\{{1_{1},1_{2},1_{3},1_{4}}\}$. The label symmetry group is just the symmetric group on four elements, whose strong generating set we can easily write down:

[TABLE]

with respect to the base $B=\langle{1,2,3,4}\rangle$. Again, re-ordering the base can be done by conjugation, thus avoiding any more expensive group theory algorithms.

3.2 The slot symmetry group

The slot symmetry group $S$ is less well-behaved, because in principle it can be any permutation group. However, Butler’s algorithm to find the canonical double coset representative only requires changing the base order of one of the subgroups. Therefore we can leave the base of the slot symmetry group intact, and change only the base of the label symmetry group.We will adopt a further convenience that the base of the slot symmetry group will be chosen to be a complete, position-ordered base. Thus if a tensor has 6 slots, the base chosen for its slot symmetry group will always be $B=\langle{1,2,3,4,5,6}\rangle$ (or, including the two columns for the sign, $B=\langle{1,2,3,4,5,6,7,8}\rangle$). Typically this results in some redundancy; for example, if a slot symmetry group only moves points $\{{1,4,5}\}$, then there is no reason to include points $\{{2,3}\}$ in the base. There is some small overhead cost in skipping past the redundant base entries; however, the tradeoff here is simplicity in describing the algorithm. We will always canonicalize with respect to strict lexicographic order $\prec$, as it is more human-readable than alternatives.121212We note that this convention is a departure from Portugal’s, which uses an ordering $\prec_{B}$ defined by the base order. So if the base is given as $B=\langle{1,4,5}\rangle$, one first extends this to a complete ordering by appending the remaining slot numbers, $B=\langle{1,4,5,2,3,6}\rangle$, and then $\prec_{B}$ is defined relative to this ordering (and hence $5\prec_{B}3$, for example). The running time of an algorithm can be sensitive to the length of the base, which in turn is influenced by the choice of base order, and in some cases it may be worth finding a more optimal choice (however, optimizing the base order is itself a costly operation!). We prefer to keep the positional/lexicographic order for simplicity, although it is not hard to modify the algorithms here for base-ordering.

3.3 The algorithm

We are now prepared to describe the Butler-Portugal algorithm for tensor canonicalization. The basic strategy is straightforward: Proceed slot by slot from 1 to $n$ and find, for each slot, the least (lexicographically) label which can be moved into that slot via a combination of slot and label symmetries. Once a slot has been filled with the least available label, that slot is frozen and we move on to the next slot. Since the slots are consumed in position-order and the base for the slot symmetry group is also in position-order, freezing a slot is the same as moving to the next $S^{(i)}$ in the stabilizer chain of $S$. The labels, on the other hand, are consumed in an order which cannot be predicted before running the algorithm. Thus we must re-order the base of the label symmetry group at every step. At each step the least possible label is placed at the first available position, and thus at the end we are guaranteed to have the arrangement of labels, out of all elements in $LgS$, that is lexicographically first.The main complication of the algorithm is that at each step, there may be multiple different ways in which the least label can be moved into the current slot, thus bifurcating the search tree. We are forced to track each possibility, because only at some future step will we be able to deduce which branch of the tree results in the lowest overall configuration of labels. For example, suppose that some tensor has a cyclic symmetry,

[TABLE]

and we wish to canonicalize an expression like

[TABLE]

where $U$ is some other tensor that has no symmetries. We can find three different ways to put the label $a$ into the first slot; one for each of the cyclic configurations of (3.7), while simultaneously exchanging the labels $(a,c)$ or $(a,b)$ as appropriate,

[TABLE]

and in the last case we have also applied the metric tensor to swap ${}^{a}{}_{a}\mapsto{}_{a}{}^{a}$. Since we must now freeze the first slot, it is clear that only the last configuration in (3.9) will lead to the least possible arrangement by lexicographic order; however, one cannot see which of these three possibilities should be kept until we examine the second slot. It is not hard to imagine situations that frustrate the decision process for arbitrarily many steps before finally allowing us to resolve the correct choice, and thus we must potentially store many intermediate configurations as we traverse the search tree. Thus in addition to iterating over the slots, we must also iterate over the set of intermediate configurations $\{{g_{i}}\}$. We store the set of configurations at the current level of the search tree in a variable Configs, and prepare the set of configurations at the next level of the search tree in Next-Configs.In Algorithm 3.1 we present a high-level view of the algorithm of [10], glossing over some details of implementation. The groups $S$ and $L$ are presumed to be stored via their bases and strong generating sets, and the base for $S$ is ordered by slot position. The algorithm will internally modify $L$ by re-ordering its base in line 19 once for each slot, when a label is selected to be placed in that slot. Each step of the main iteration over slots is split into two phases: First, in lines 4–18, we obtain the least available label that can be moved into the current slot under the action of the symmetry groups $S^{(i)}$ and $L^{(i)}$, which respectively stabilize the already-consumed slots and the already-consumed labels. At the same time, we populate the set Config-Slots-of-Least-Label with ordered pairs $(g,j)$ where $g\in\textsf{{{Configs}}}$ is one of the nodes at this level of the search tree, and $j$ is a slot reachable from $i$ (via some $s\in S^{(i)}$) wherein it is possible to encounter the least-available label (via some $\ell\in L^{(i)}$). In the second phase, from lines 19–25, we iterate through each of the $(g,j)$ just collected and populate the set Next-Configs with the next level of the search tree, given by $\ell\llbracket g\llbracket s\rrbracket\rrbracket$, where $s\llbracket i\rrbracket=j$ and $\ell\llbracket g\llbracket j\rrbracket\rrbracket=\textsf{{{global-least-label}}}$.The function \textsc{Re-order-Base}\big{(}L,i,\textsf{{{global-least-label}}}\big{)} in line 19 re-orders the base of $L$ such that the $i$-th base point becomes global-least-label, while keeping the first $i-1$ points unmodified. The main insight of [10] is that Re-order-Base can be implemented for the label symmetry group via conjugation by the appropriate label swap. The function \textsc{Coset-Rep}\big{(}j,S^{(i)}\big{)} in line 22 returns a group element of $S^{(i)}$ which maps slot $i$ to slot $j$ (there may be several such group elements; it does not matter which is returned). Similarly, the call to \textsc{Coset-Rep}\big{(}g\llbracket j\rrbracket,L^{(i)}\big{)} in line 23 returns a group element of $L^{(i)}$ which maps the $i$-th base point (which was earlier set to global-least-label) onto the label $g\llbracket j\rrbracket$. The functions Permutation-Inverse and Append are straightforward.We now comment briefly on the data structures needed to represent the various quantities in Algorithm 3.1. In order to implement Coset-Rep we presume that some sort of Schreier tree structure (various examples are described in [16, 17, 18, 19]) has been constructed in the course of enumerating the orbits $\Delta^{S}_{i}$ and $\Delta^{L}_{g\llbracket j\rrbracket}$; Coset-Rep then constructs a group element by tracing the Schreier tree from the specified leaf back to the root. In fact, since the base of $S$ is never re-ordered, it is possible to compute $\Delta^{S}_{i}$ ahead-of-time, if computational time is at a higher premium than memory space. There are a few different options to store all of the $\Delta^{S}_{i}$ together with their Schreier trees, such as a labelled branching [18] or a Schreier vector structure [19].131313In our particular implementation of the algorithm in Section 4.3, we have chosen a “cube Schreier tree” structure as detailed in [17]. This structure requires more time to create than a standard Schreier tree, but the tree is guaranteed to be balanced and thus minimizes the time used by Coset-Rep. There are also other small optimizations, such as for the intermediate array Config-Slots-of-Least-Label to store not the actual configuration elements $g$, but rather an index into the Configs array, etc.

3.4 Worst-case complexity analysis

In [10], Portugal shows by direct testing that for the case of many Riemann tensors with random contractions, Algorithm 3.1 appears to have an average-case time complexity of $O(n^{5})$ for $n$ the number of slots (in fact, our measurements in Section 5 show that the particular implementation in Mathematica runs slightly faster then $O(n^{4})$). While the Riemann tensor is presented as having one of the more complicated slot symmetry groups that is typically encountered in computations, in fact it represents one of the easier (non-trivial) cases for the algorithm to resolve. The worst-case behavior actually comes from fully symmetric tensors (or fully antisymmetric) when they are fully contracted, for example:

[TABLE]

Such tensor monomials suffer from a re-labelling ambiguity which maximally frustrates the decision-making of the algorithm, rendering it completely unable to prune the search tree. Let us follow this example through a few iterations of the algorithm to see why:First, we examine slot 1. There are 6 possible ways to put the least label $a$ into slot 1, which correspond to the following slot and label symmetries:

[TABLE]

which in turn give the configurations

[TABLE]

Now we must examine slot 2 on each of these configurations. For each configuration, we will find that there are 5 distinct ways to move label $b$ into slot 2, and thus the set of configurations will grow to $6\cdot 5=30$ entries. At the next step, on each of these entries we will find 4 ways to move label $c$ into slot 3, and so on, so that the largest intermediate result contains $6!=720$ entries, and the total time examining all of the intermediate results is $1+6+6\cdot 5+6\cdot 5\cdot 4+\ldots+6!=1957$, multiplied of course by whatever time it takes to execute each iteration of the internal loops. Clearly in such cases the Butler-Portugal algorithm will experience combinatorial explosion in both time and space. In a sense this behavior is especially frustrating, because to the end-user with human eyes, the canonical result is “obviously” the one with all of the indices sorted:

[TABLE]

and yet the algorithm will take quite a while to work this out. Increase the number of indices to 12 and the algorithm will use 24 GB of RAM and take a full day. Increase the number of indices to 14 and the algorithm (if allowed to write intermediate results to disk) will completely fill a 4 TB hard drive over the course of 30 years.141414The space requirements here are calculated assuming each configuration is an array of $(n+2)$ 32-bit integers. The timings are projected from a best-fit curve of the test in Figure 8(a) involving up to 10 indices, which thankfully run in only a minute or two per test iteration (on a machine with a 2.40 GHz processor and 8 GB RAM). Note that it is difficult to get a reliable fit to these curves due to the paucity of data points, and thus the projected timings may be wildly inaccurate.This problem is noted by the authors of [5], and they attempt to mitigate it by first sorting any tensor products (such as $TS$ here) by increasing group order of their slot symmetry groups. Thus $S$, since it has no symmetries, will be put first, and then the algorithm becomes linear as there is no longer any ambiguity of intermediate results. This is a dramatic improvement; however, it is still thwarted when both tensor factors have a high degree of symmetry, as in

[TABLE]

Products of tensors with many symmetric or antisymmetric slots show up in several contexts, such as in higher-spin theory [20, 21], supergravity [22, 23], representation theory [24, 25, 26], and the construction of spherical harmonics with high quantum numbers and/or in higher dimensions based on homogeneous polynomials [27]. Therefore it is worth resolving this problem more thoroughly.

4 The improved algorithm

If one is doing the sort of computation where one expects to see a lot of symmetric/antisymmetric tensors, even of with a moderate number of indices, the behavior of Section 3.4 is troubling. One needs an efficient solution. One naïve strategy is to treat fully (anti)symmetric tensors in a special way, perhaps by simply sorting their indices (and inserting the appropriate sign for antisymmetric tensors). But we can quickly see that this will not be sufficient, for two reasons: First, such tensors can appear in tensor products with tensors with other symmetries, and thus to be fully general, an algorithm must be able to deal with arbitrary (anti)symmetric subsets of indices, rather than having special code only for the case of full (anti)symmetry. Second, when including dummy contractions in such tensor products, it is not enough merely to sort the (anti)symmetric subsets of indices, because dummy labels may be exchanged while canonicalizing neighboring tensors, thus causing the notion of “sorted” to change as well. A more intelligent approach is needed.

4.1 Insights from Penrose graphical notation

It turns out that a useful way to approach the problem is to employ the Penrose graphical notation for tensor contractions, first published in [28] (see also [29] for an alternative version). In this notation, tensors are represented by various (arbitrary) shapes, and their indices are represented by lines emanating from these shapes. For example, one might have tensors $\theta$ and $\chi$ given by:

[TABLE]

Lines going to the top of the diagram represent upper indices, while lines going to the bottom of the diagram represent lower indices. To represent contracted indices, we merely connect the lines:

[TABLE]

Of course, the label $f$ is meaningless; the important information is the line connecting $\theta$ and $\chi$ which shows the contraction. One also has special symbols for the metric, inverse metric, and Kronecker delta,

[TABLE]

which are made only of connecting lines, and which obey the obvious graphical relation

[TABLE]

Combining these with the shapes for tensors like $\theta$ and $\chi$, one can create diagrams for tensor contractions of arbitrary complexity. There are some additional features of the notation which we will not need here; we refer the reader to [28] for more detail.

4.1.1 Dummy pairs as links for symmetry propagation

For our purposes, it is useful to enlarge the diagrams a bit and write in the slot number to which a line is attached. For example, the tensor contraction

[TABLE]

can be represented by the diagram in Figure 1, where the slots are numbered 1-12 in the order they appear in (4.5).

It is instructive to start from the arrangement in Figure 1 and walk through the steps a human might take to canonicalize it. Suppose, for example, that the tensor $T$ is totally symmetric over the slots $\{{3,4,5,6}\}$. This situation is depicted in Figure 2(a). The dummy labels $c,d,e,f$ are not meaningful, but the objective is to find a combination of slot- and label-rearrangements which minimizes (4.5) with respect to lexicographical ordering. The first step is to “disentangle” the contracted edges in Figure 2(a) by moving the endpoints which are attached to the slots $\{{3,4,5,6}\}$ where the slot symmetry acts. This results in

[TABLE]

as shown in Figure 2(b). The final step is to re-name all of the dummy labels in order of appearance, which produces

[TABLE]

as shown in Figure 2(c). Assuming the tensor $U$ has no slot symmetries, then we have achieved the least possible lexical order, with the dummies $\langle{c,d,e,f,c,d,e,f}\rangle$ repeated in order.The above process is a strict application of the available symmetry groups: first using slot symmetries to make the two sets of dummies match in order $\langle{e,d,f,c,e,d,f,c}\rangle$, and then using label symmetries to rename these into the correct order. It is easy to describe this process in words, but to put it in an algorithm which visits each slot only once, scanning from left to right, is a bit non-trivial. Since the Butler-Portugal algorithm does not backtrack once it has made a decision to place a label in a slot, it must instead retain enough information that it can make the decision later. In this case, it will see $4!$ ways to put $\langle{c,d,e,f}\rangle$ into slots $\langle{3,4,5,6}\rangle$, and only begin to resolve which of those possibilities to retain when it reaches slot 7. The decisions needed to reduce $4!$ possibilities down to the 1 correct one are not fully made until it visits slot 10.

Of course, we do not wish to introduce backtracking into the algorithm. Instead, we will slightly redefine the problem in a way that makes it amenable to a strict left-to-right pass without having to retain (and iterate over!) a factorial number of intermediate states. Our main observation is that the symmetries at one end of an edge in a Penrose tensor diagram induce symmetries of the opposite end. Thus the slot symmetries of tensor $T$ can “propagate” along the connecting lines to tensor $U$, and induce new symmetries which appear to act on the slots of $U$ (as we shall see, however, they actually act on the labels!). This fact is somewhat obvious from the Penrose graphical notation, but is obscured in the standard index notation because the “propagation” of symmetries along contractions may be effected only by executing slot and label symmetries at the same time.To illustrate, we walk through the same example. We begin with (4.5), as shown in Figure 3(a). Now, instead of rearranging the endpoints attached to $\{{3,4,5,6}\}$, let us actually freeze them in place; after all, the labels there are already $\langle{c,d,e,f}\rangle$, which is what we want. But there is a symmetry among the slots $\{{3,4,5,6}\}$ which we have not used, and which is crucial to obtain the correct result. Let us “propagate” this symmetry forward, along the legs attached at $\{{3,4,5,6}\}$, to the slots $\{{7,8,9,10}\}$ on tensor $U$. In fact, for reasons we will soon explain, we will consider this symmetry to be associated to the labels $\{{e,d,f,c}\}$ which are present at slots $\{{7,8,9,10}\}$, as depicted in Figure 3(b).

The tensor $U$ does not actually have a symmetry in the slots $\{{7,8,9,10}\}$, but by simultaneously applying the slot symmetries on $\{{3,4,5,6}\}$ and the standard label symmetries, one can give the appearance that $U$ has such a symmetry:

[TABLE]

This is what we mean by a “propagated” symmetry. Thus we may freeze the $T$-attached endpoints of the contractions in place, and treat the $U$-attached endpoints as though they are free to move. This allows us to rearrange the labels $\{{e,d,f,c}\}$ on tensor $U$ independently of those on tensor $T$, and thus reach the canonical configuration as shown in Figure 3(c).With this notion of symmetry propagation along dummy contractions, we can effectively implement the necessary decision process in a single left-to-right pass without storing many intermediate configurations. Instead we must store information about the symmetries to be propagated forward. For arbitrary symmetry groups, this might become unwieldy; but if this method is used only for totally (anti)symmetric subgroups, then it is simple to implement. One needs only a single array of length $n$ to store all of the information regarding propagated symmetries (as will be described in Section 4.2). One can then take advantage of this information to make a single definite decision at each new visited slot, completely eliminating the $O(n!)$ intermediate steps needed in the Butler-Portugal approach.We have stated that slot symmetries which have been propagated along dummy contractions become label symmetries, and this requires some explanation. Up until now, the tensor $U$ had no symmetries of its own. But suppose instead it had a single symmetry $(9,11)(10,12)$ which exchanges the last two pairs of slots. Our lexicographical ordering prioritizes free indices ahead of dummies, so we should use this symmetry to move $\langle{g,h}\rangle$ ahead of $\langle{f,c}\rangle$. But the propagated symmetries, which originated from total symmetry of slots $\{{3,4,5,6}\}$ of $T$, must remain associated to the contraction edges which attach to slots $\{{3,4,5,6}\}$. Therefore, the propagated symmetries will now be associated with slots $\{{7,8,11,12}\}$ as shown in Figure 4. It is clear that we should really associate the propagated symmetry with the labels $\{{e,d,f,c}\}$, and thus the generators

[TABLE]

should be appended to the label group $L$ (here, as before, the subscript 2 refers to the appearance of the label as an upper index; namely on the tensor $U$).

4.2 Data structures

In this section we give several data structures which are needed to implement the symmetry-propagation mechanism described in Section 4.1.1. These data structures are also useful for improving the basic Butler-Portugal algorithm as they can reduce most of the label-group computations of Algorithm 3.1 to mere array lookups. These data structures we will refer to as the Values array, the Label-Groups array, the Symmetric-Subsets array, and the Propagated-Symmetries array. All of these are simple length-$n$ arrays, and thus using these data structures incurs an additional $O(n)$ memory cost only.We will also describe Config-Slots-of-Least-Value, which is a slightly different way of tracking the current level of the search tree, and the notion of a Least-Value-Set, which is just a useful way to group the entries of Config-Slots-of-Least-Value.

4.2.1 The Values array

Lines 4–18 of Algorithm 3.1 are occupied with finding the least available label which can be moved into the current slot. The costliest parts of this code are line 4 and line 8, where the orbits $\Delta_{i}^{S}$ and $\Delta^{L}_{g\llbracket j\rrbracket}$ are computed. In fact, the slot orbit $\Delta_{i}^{S}$ can be computed once-and-for-all at the time a tensor and its symmetries are declared, and stored in a type of Schreier tree structure, which we will not detail here (see, for example, [16, 17, 18, 19]). However, the label orbit $\Delta^{L}_{g\llbracket j\rrbracket}$ depends on which labels have been placed in a tensor’s slots, which comes from the context of a specific tensor expression. It would appear that this must be computed on-the-fly.However, we can introduce the concept of the value of a label, which can achieve two useful things simultaneously: First, it can reduce lines 8 and 9 of Algorithm 3.1 to a simple array lookup, thus skipping the computation of $\Delta^{L}_{g\llbracket j\rrbracket}$ altogether; second, it can concisely encode all the data as to which indices are free, or belong to a given dummy set or repeated set, into a single array of length $n$. The concept is simple, and satisfies the properties of Definition 4.2.1:{definition}Properties of the Values array.

1. V1.

The value of a label is equal to the (number encoding) the least label in the group to which a label belongs; 2. V2.

Each label in a group of interchangeable labels has the same value; 3. V3.

Every free index belongs to its own group; 4. V4.

A group of interchangeable labels is usually contiguous in label space.151515The exception is when a label belongs to a dummy set for a vector bundle on which there is no defined metric. Then upstairs labels cannot mix with downstairs ones. The mechanics of this case will be explained in Section 4.2.2.

Take, for example, the tensor expression

[TABLE]

where as before, we insert the subscripts ${}_{1},\,_{2}$ to distinguish identical labels. The labels have a natural lexicographical ordering $\langle{a,c,1_{1},1_{2},b_{1},b_{2}}\rangle$ which we use as the index into the array. Then the Values array can be written

[TABLE]

As a more elaborate example, consider again the contraction

[TABLE]

The Values array corresponding to this expression is given by

[TABLE]

Notice in each case that all labels which can be interchanged by label symmetries have the same value; moreover, that value is equal to the (numerical) index of the least label itself. Then, for example, if one needs to know the least label that can be moved into slot 7 of (4.12), the process is simple:

Look up the label in slot 7 of (4.12): $e_{2}$. 2. 2.

In the Values array (4.13), look up the value of $e_{2}$: 5. 3. 3.

In the list of sorted labels $\langle{a,b,g,h,c_{1},c_{2},d_{1},d_{2},e_{1},e_{2},f_{1},f_{2}}\rangle$, take the 5th entry: $c_{1}$.

and thus we see that $c_{1}$ is the least label that can be moved into slot 7, without having to explicitly enumerate the orbit $\Delta^{L}_{g\llbracket 7\rrbracket}$.Of course, since the algorithm visits the slots in order, we should have already placed $a$, $b$, and four of the dummy labels by the time we reach slot 7. Once a label has been consumed by placing it in a slot, it is no longer interchangeable with the unused labels, and the Values array must reflect that fact for the correct functioning of the algorithm. As we shall detail in Section 4.3, the Values array will be updated at the end of every step to ensure that the properties of Definition 4.2.1 continue to hold.

4.2.2 The Label-Groups array

The Values array can efficiently keep track of which subsets of labels are interchangeable; however, it does not contain information about how they can be interchanged. In Algorithm 3.1, this information is contained in the label symmetry group $L$, which must be continually updated in line 19 via the procedure Re-order-Base whenever any two labels are swapped. The main insight of [10] is that Re-order-Base can be implemented by conjugating the members of $L$, which is much faster than a generic base-change algorithm.However, even this can be improved upon. To this end, we introduce another data structure of length $n$, the Label-Groups array. Instead of literally storing a base-and-strong-generating-set for the label group $L$, we will store an array of group codes which indicate how labels may be exchanged among a given equal-value set. Like the Values array, the Label-Groups array is indexed by the labels, in lexicographic order, and thus should be thought of as living “in label space”. Each entry in the Label-Groups array shall be one of the group codes as described in Table 1 (these can be implemented, for example, by an enum type or similar).

It is useful to display the Values array along with the Label-Groups array, as together these constitute all the information about the label symmetry group $L$. To revisit our earlier examples, the expression

[TABLE]

corresponds to the Values and Label-Groups

[TABLE]

and the expression of our earlier example,

[TABLE]

corresponds to

[TABLE]

if we assume that there is a metric which can exchange upper and lower indices. If no metric exists, then we must remember that only labels which can be exchanged can have the same value; thus the Values and Label-Groups arrays will read

[TABLE]

In this case, the set of labels with value 5 is not contiguous (likewise the set with value 6), and thus our algorithm cannot strictly rely on the Values array being a non-decreasing sequence. However, the entries in the Label-Groups array give enough hinting to deduce that the Values array must have a contiguous sequence of $\langle{5,6}\rangle$ blocks, and this fact we can take advantage of.In any of these cases, the Values and Label-Groups arrays in (4.15), (4.17), and (4.18) contain sufficient information to reconstruct any necessary label exchange to bring the least label into a slot containing a given label. Suppose, for example, we are given the label $e_{2}$ in the presence of a symmetric metric, thus corresponding to (4.17). We construct the label exchange as follows:

The value of $e_{2}$ in (4.17) is 5, which means it can be swapped with the label in the 5th position, $c_{1}$. 2. 2.

The group of $e_{2}$ is $\mathord{+}$, which means that each of the swaps $(c_{1},e_{1})(c_{2},e_{2})$ and $(e_{1},e_{2})$ are allowed. 3. 3.

The position of $e_{2}$ is 10; since $(10-5)$ is odd, we determine that both exchanges must be made in order to map $e_{2}\mapsto c_{1}$. This can be achieved by

[TABLE]

By contrast, suppose instead that we have no metric, as in (4.18). In this case, the label exchange can be constructed as follows:

The value of $e_{2}$ in (4.18) is 6, which means it can be swapped with the label in the 6th position, $c_{2}$. 2. 2.

The group of $e_{2}$ is $\mathord{\top}$, which means it is an upper index. In our conventions, the lower indices are ordered before each corresponding upper one; therefore we conclude that each of $c_{2},\,e_{2}$ should be paired with the label immediately to its left in the array. 3. 3.

With this information, we can construct the necessary label exchange: $(c_{1},e_{1})(c_{2},e_{2})$.

Thus we can see how the Values and Label-Groups arrays can be used to quickly circumvent many of the more expensive steps in the inner loops of Algorithm 3.1. It should also be straightforward how one might construct these arrays in the first place; the Values are merely a means of encoding the various sets of indices that appear in an expression, and the Label-Groups are fixed entirely by the combination of label types and type of metric tensor currently defined (or not defined) on a given vector bundle.

4.2.3 The Symmetric-Subsets array

The previous data structures can be used to improve the performance of the standard Butler-Portugal algorithm even without considering the new concepts introduced in Section 4.1.1. In this section and the next, we will introduce further data structures for implementing the symmetry-propagation ideas which allow fast canonicalization of totally (anti)symmetric subsets of indices.The first of these is the Symmetric-Subsets array, which is an extra input to the algorithm, of length $n$, which stores condensed information on where to find (anti)symmetric subsets in the slot symmetry group $S$. The Symmetric-Subsets array can be generated by $O(n^{2})$ group-membership tests (slightly modified for signed permutations) to find the elementary length-2 cycles which generate each (anti)symmetric subgroup (see [9, 13, 14]). Since the information contained in Symmetric-Subsets is derived only from the slot symmetry group $S$, it can be obtained ahead of time at tensor declaration time and stored at $O(n)$ cost, which is trivial compared to the cost of the Schreier vector structure needed for $S$ itself.161616We point out that the Symmetric-Subsets array in fact describes redundant information which is already strictly present in the Schreier vector structure of $S$. However, this information is in a special format which allows quick manipulation by taking advantage of the simple structure of totally (anti)symmetric groups. One could set Symmetric-Subsets to the zero array and the algorithm will still produce correct results, but more slowly.{definition}Properties of the Symmetric-Subsets array.

1. S1.

A zero entry indicates that a slot does not participate in any (anti)symmetric subsets. 2. S2.

Positive integers represent a symmetric subset of slots. All slots with the same positive integer can be exchanged with each other. 3. S3.

Negative integers represent antisymmetric subsests. All slots with the same negative integer can be exchanged at the cost of a minus sign. 4. S4.

In addition, we shall require that the integers be chosen in a sequence of increasing absolute value, and that each absolute value is used for only one subset (thus, 1 and $-1$ should not both appear, but rather 1 and $-2$, etc.).

The structure of the Symmetric-Subsets array is simple: The positions in the array represent slots (and thus Symmetric-Subsets lives in slot space, in contrast to Values and Label-Groups which live in label space). Each entry is an integer, with meanings as given in Definition 4.2.3. For example, suppose we have a tensor $T_{abcdef}$ which is symmetric on the last 4 slots; thus the slot symmetry group $S$ is given by the strong generating set

[TABLE]

The corresponding Symmetric-Subsets array is then given by

[TABLE]

Similarly, the Riemann tensor $R_{abcd}$ has the slot symmetry group given by the strong generating set

[TABLE]

and the corresponding Symmetric-Subsets array is

[TABLE]

Finally, consider the contraction $T_{abcdef}R^{cd}{}_{gh}$. This has the slot symmetry group given by

[TABLE]

and the Symmetric-Subsets array will read

[TABLE]

Of course, we know that the expression $T_{abcdef}R^{cd}{}_{gh}=0$ because the contraction of the indices $cd$ will bring the symmetry of slots $\langle{3,4}\rangle$ in conflict with the antisymmetry of slots $\langle{7,8}\rangle$. Thus we can appreciate the usefulness of having pre-computed the structure in (4.25); in cases like this one, this will allow our algorithm to detect zero contractions much earlier than Butler-Portugal, as will be shown in Section 5.In addition we note that Symmetric-Subsets, being a container for slot-symmetry information, does not take into account that $cd$ are dummy indices while $ef$ are free; the slots $\{{3,4,5,6}\}$ participate in the same slot symmetry regardless of which indices they contain.

4.2.4 The Propagated-Symmetries array

In this section we will discuss the Propagated-Symmetries array, whose purpose is to track the information needed to implement the symmetry-propagation mechanism of Section 4.1.1. The types of information that need to be stored are twofold: First, one needs to know which portions of the Symmetric-Subsets are attached to dummy indices; second, one needs to know where the other “ends” of those dummy contractions reside, in the sense of the Penrose graphical notation introduced in Section 4.1. Both of these needs can be served using a single array of length $n$, whose structure is similar to that of Symmetric-Subsets, but with an additional distinction between odd- and even-numbered entries, which will be explained shortly.Before defining the entries in Propagated-Symmetries, however, we must decide whether it should be indexed by slots or labels. We first note that the symmetries described in Propagated-Symmetries are label symmetries, as suggested by Figure 4, and thus it seems sensible for Propagated-Symmetries to live in label space. However, as will be explained in Section 4.2.4, the Propagated-Symmetries array can never be accessed directly, because each configuration $g\in\textsf{{{Configs}}}$ which we must iterate over (as in line 6 of Algorithm 3.1) may have been reached by a different combination of slot and label symmetries. That is, the labels may have been renamed, and/or the slots shifted. The total configuration is given by $g=\ell\llbracket g_{\text{init}}\llbracket s\rrbracket\rrbracket$ where $(\ell,s)$ are the independent label- and slot-actions which have been accumulated in the course of the algorithm.171717In our version of the Butler-Portugal algorithm in Algorithm 3.1, we store only the total configuration $g$, whereas in the original version in [10], the label- and slot-actions $(\ell,s)$ are stored separately. The Propagated-Symmetries array is a global variable whose information is needed in every configuration we iterate over. Thus if Propagated-Symmetries is label-indexed, we must know $\ell$ in order to access it properly; alternatively, if it is slot-indexed, we must know $s$. Neither convention seems to offer a computational advantage. We will choose Propagated-Symmetries to be slot-indexed, despite that fact that its entries represent label symmetries, because we find it easier to relate to the diagrams of Section 4.1 this way.{definition}Properties of the Propagated-Symmetries array.

1. P1.

A zero entry indicates that the label at a given slot does not participate in any (anti)symmetric subsets. 2. P2.

An odd integer (positive or negative) indicates that the label at a given slot has a symmetry or antisymmetry, according to the sign of the integer. All labels at slots with the same integer participate in the same symmetry. Odd integers may be entered for both dummy labels and component labels. 3. P3.

An even integer (positive or negative) indicates that the dummy label at a given slot has a symmetry/antisymmetry which has been propagated from its corresponding partner with odd-numbered entry. If the even-numbered entry is $k$, then its partner symmetry is the odd number of the same sign, and one less in absolute value, $(\operatorname{sign}k)\cdot({\lvert{k}\rvert}-1)$. 4. P4.

The process which creates the entries shall work left-to-right in slot order, alternatively enumerating the odd-numbered (anti)symmetric sets and then propagating to their corresponding even-numbered partners. If at any time this process would cause a new (nonzero) entry to overwrite an old (nonzero) entry, the entry with lower absolute value shall take precedence (such conflicts can happen when both ends of a contraction are attached to the same set of symmetric slots, or when the “far” end of a contraction is attached to another symmetric subset, distinct from the subset to which the “near” end is attached). 5. P5.

If there are labels whose slots have the same Symmetric-Subsets entry, but whose Values entries are different (thus, they belong to different component-label or dummy sets), then they will be given different odd integers; likewise, any of their propagated partners must receive different (and corresponding) even integers. Thus the entries in Propagated-Symmetries must always represent valid label symmetries. 6. P6.

Nonzero entries of either type will be entered only if they are shared by at least two dummy labels (after applying Rule P4); “symmetric subsets of length 1” are equivalent to non-symmetric labels, and should have the entry 0. 7. P7.

As with the Symmetric-Subsets array, odd integers will be chosen in a sequence of increasing absolute value, and each absolute value will be used for only one subset. Even integers are always entered as partners of odd ones, and hence also form such a sequence. The even-numbered sequence may have gaps, either because component labels do not have partners, or as a result of Rule P4.

The properties satisfied by Propagated-Symmetries are given in Definition 4.2.4. Of all of the data structures we will introduce, this one has the most complex system of rules, and some examples will help clarify how it is to be filled out. First, suppose we have a tensor $T_{abcdef}$ which is symmetric in its last 4 slots, and thus has the Symmetric-Subsets array

[TABLE]

We have previously considered the expression

[TABLE]

whose Values array is

[TABLE]

Now consider how to fill out its Propagated-Symmetries array according to the rules in Definition 4.2.4. It has one subset of equivalent component labels $\{{1_{1},1_{2}}\}$, but they are in slots $\langle{1,2}\rangle$ which do not have a slot symmetry. Therefore the $\langle{1,2}\rangle$ entries of Propagated-Symmetries should be 0. Next, there is one dummy pair, $\{{b_{1},b_{2}}\}$ in slots $\langle{4,5}\rangle$. These slots are in a symmetric subset with each other, and thus should receive the odd-number entry 1. Because of the precedence established in Rule P4, there will be no entries for propagated symmetries, because both ends of the $b_{1}\smile b_{2}$ contraction are already attached to the same symmetric subset. The Propagated-Symmetries array is therefore

[TABLE]

Next, consider a slightly more elaborate example: the expression $T_{abcdef}R^{cd}{}_{gh}$. Again we take $T$ to be symmetric on its last 4 slots, and $R$ is the Riemann tensor. We first write out the Values array (remembering that free indices are always sorted before dummies):

[TABLE]

Next, the Symmetric-Subsets array is as in (4.25):

[TABLE]

Finally, to construct the Propagated-Symmetries array, we note that the only non-free labels are in slots $\{{3,4,7,8}\}$, and by the precedence given in Rule P4, we obtain

[TABLE]

In particular, we do not get $-3,-3$ in slots $\langle{7,8}\rangle$. This we see that Propagated-Symmetries allows us to “look ahead” and see that this contraction must be zero, because the entries in slots $\langle{7,8}\rangle$ of Propagated-Symmetries disagree in sign from those in $\langle{7,8}\rangle$ of Symmetric-Subsets.As a final example, consider the expression $T_{ab11cd}R^{c}{}_{e}{}^{d}{}_{f}$, again with the same tensors. This time, the Values array is

[TABLE]

and the Symmetric-Subsets array remains the same as in (4.25):

[TABLE]

However, to construct the Propagated-Symmetries array, we must apply Rule P5, since there are both component labels and dummy labels appearing in the symmetric subset of slots $\{{3,4,5,6}\}$, and these label sets have different entries in the Values array. The symmetric subset “splits” into two pieces, giving the Propagated-Symmetries array

[TABLE]

Notice that slots $\{{7,9}\}$ must be marked $4,4$, which shows that they have an induced symmetry propagated from slots $\{{5,6}\}$. There are no entries with the number 2, because the labels in slots $\{{3,4}\}$ are component labels and do not have partners.Finally, note that because we have chosen Propagated-Symmetries to be slot-indexed, we must store the slot symmetry $s$ which was used to reach the current configuration $g$ from the initial one $g_{\text{init}}$. In Algorithm 4.1, we will re-define Configs to be the set of ordered pairs $(g,s)$. Alternatively, one could store $(\ell,s)$ as was done in the original implementation of Butler-Portugal [10] and use the relation $g=\ell\llbracket g_{\text{init}}\llbracket s\rrbracket\rrbracket$; however, we find it more useful to store $(g,s)$ as $\ell$ alone is not typically needed.181818If we had chosen Propagated-Symmetries to be label-indexed, we would instead need to store $(g,\ell)$.

4.2.5 The Config-Slots-of-Least-Value array and its Least-Value-Set entries

One final data structure we will need is the Config-Slots-of-Least-Value array, which replaces the Config-Slots-of-Least-Label array appearing in Algorithm 3.1. The overall purpose of Config-Slots-of-Least-Value will be the same: to store the configurations and slot numbers at which one can find a label with least value. However, because the properties of the label symmetry group are now distributed between three data structures Values, Label-Groups, and Propagated-Symmetries, some slight changes will be necessary.First, as mentioned in Section 4.2.4, the Configs array will now have to store ordered pairs $\{{(g,s)}\}$ where $g$ is the total slot-to-label map of the configuration, and $s$ is the total action of the slot symmetry group which brought us to $g=\ell\llbracket g_{\text{init}}\llbracket s\rrbracket\rrbracket$. Therefore each entry in the Config-Slots-of-Least-Value array must contain a reference to the ordered pair $(g,s)\in\textsf{{{Configs}}}$ from which the current least-value slot numbers originate (that is, $(g,s)$ is a reference to the parent node in the search tree).In addition, each entry of Config-Slots-of-Least-Value must contain an ordered pair slot numbers $(p,q)$ rather than just a single one. The reason for this is because label symmetries can come from two different sources: the original label symmetry group $L$ which is represented by the combination of Values and Label-Groups; and the new label symmetries which come from propagation of slot symmetries, as represented in Propagated-Symmetries.191919One could alternatively use a Schreier-tree structure to store the entire label group as was done in Algorithm 3.1, but this would sacrifice quite a bit of speed, as one would have to add the propagated symmetries by “sifting” them into the Schreier tree, and then re-compute orbits, etc.Thus the Config-Slots-of-Least-Value array will have the following structure:

[TABLE]

where one can think of each $(g_{k},s_{k})$ as a reference into the Configs array. We do not expect each of the $(g_{k},s_{k})$ to be distinct, as it is possible that a given configuration contain multiple labels of (equal) least value in the orbit of the slot currently under consideration. We will find it useful to collect together all of the Config-Slots-of-Least-Value coming from the same configuration $(g,s)$ into a Least-Value-Set; thus, the Least-Value-Set corresponding to a given $(g,s)$ is precisely the set of child nodes of the search tree under the parent node $(g,s)$.202020Thus Config-Slots-of-Least-Value could be implemented as a multi-dimensional array containing Least-Value-Set’s, although one must take into account that each Least-Value-Set may be of different length. For the timing data in Section 5, we have chosen a flat implementation. The meanings of the slot pairs $(p_{k},q_{k})$ are as follows:

The $p_{k}$ are the slot numbers in the slot-symmetry orbit $\Delta^{S}_{i}$ where a least-value label can be found directly. 2. 2.

The $q_{k}$ are the slot numbers where a least-value label can be found by application of Propagated-Symmetries, even if that slot is outside the orbit $\Delta^{S}_{i}$.

Usually, one has $p_{k}=q_{k}$ for each $k$, except in the case where the Propagated-Symmetries have been used to step outside the usual slot-symmetry orbit to find the least-value label. Then $p_{k}$ gives the point in the slot-symmetry orbit we jumped away from, and $q_{k}$ gives the place where the least-value label was actually found. These two slot numbers together give us implicit information about how to reconstruct the necessary slot- and label-symmetries to put the least label into our current position.

4.3 The main algorithm

Having developed the necessary preliminary concepts, we now present the improved algorithm. We will first give an overview in Algorithm 4.1 of the main loop which iterates once over the slots of the initial configuration $g_{\text{init}}$. There are several subprocedures called by the main loop, most of whose code we relegate to Appendix B, with the exception of Append-Non-Redundant-Instances which contains the essential logic that prevents iteration over redundant branches of the search tree.

The general strategy of Algorithm 4.1 is the same as Algorithm 3.1: namely, to visit each slot one at a time from left to right, determine the least label which can be moved into that slot, place that label, and then move on to the next slot. As before, if there are multiple ways in which the (same) least label can be moved into the current slot, then the search tree bifurcates, and the next iteration must look at all of the resulting possibilities. The key difference is that in Algorithm 4.1, we make an effort to prune these extra branches whenever they would be redundant. Not all types of redundancy are detected, but we do handle the most common case of (anti)symmetric subgroups filled with dummy labels. We apply the symmetry propagation concepts of Section 4.1.1 in order to make an early decision about placing such dummy labels into slots, while preserving the symmetry information needed to ensure that later slots can be canonicalized. The same information can be used for early detection of zero results.The main loop over the slots of $g_{\text{init}}$ runs from lines 4 through 21 of Algorithm 4.1. In line 5, we obtain the orbit $\Delta^{S}_{i}$ of the current slot $i$ under the slot symmetry group $S$ (an operation which can be made quite fast if we have stored $S$ in the form of a Schreier tree structure ahead of time). In line 6, we populate the Config-Slots-of-Least-Value array described in Section 4.2.5, as well as record the least-value label which can be moved into this slot (either via the slot symmetry group $S$, or via the extra label symmetries in Propagated-Symmetries).Next in lines 8 through 15, we iterate over the next level of nodes of the search tree, grouped by their parent node $(g,s)\in\textsf{{{Configs}}}$ in the form of a Least-Value-Set. On each Least-Value-Set we make two passes; first to iteratively add information to the Propagated-Symmetries array, and then to use the information to make decisions about which leaf nodes to append to Next-Configs. Note that the data in Propagated-Symmetries is filled out according to the rules in Definition 4.2.4, but it is not done all at once; rather, only the Symmetric-Subsets which are encountered within the current Least-Value-Set are entered. Thus at any given time before the last iteration of the algorithm, the Propagated-Symmetries array is incomplete; however, it always contains just enough information to be used in Append-Non-Redundant-Instances.212121Propagated symmetries always come from slot symmetries, and if there is a symmetric subgroup to be propagated, then it is always true that all of its slot numbers should be listed in the orbit $\Delta^{S}_{i}$, and thus appear in the current Least-Value-Set. Thus at the very latest, each symmetry which can be propagated is entered into Propagated-Symmetries immediately before it is needed for the decision-making phase. One could populate the entire Propagated-Symmetries ahead of time, but doing so iteratively as we do here prevents us having to do any more work than necessary (in the event, for example, that the entire algorithm is short-circuited by a zero detection).In lines 16 through 20, we do some final steps before moving on to the next slot iteration. First we must update the Values and Label-Groups arrays to maintain the properties in Definition 4.2.1 and Table 1. This effectively implements “re-ordering the base” for the label symmetry group, as the lowest-value label will be marked as used and its exchange symmetry with other labels will be erased. Finally we sort and remove duplicates from Next-Configs, copy the result back into Configs, and check for zero.The subprocedures Get-Least-Value-Instances, Update-Propagated-Symmetries, and Update-Values-and-Label-Groups are intended to populate the various data structures defined in Section 4.2 and maintain the conditions of their definitions for each loop iteration. The conditional Zero-Due-to-Propagated-Symmetries constitutes the early check for zero using the extra information assembled in the data structures of Section 4.2. However, we emphasize here that this check is not optional; rather, it is necessary to take into account all the possibilities for a zero result that might be contained in Propagated-Symmetries immediately. This allows Append-Non-Redundant-Instances to be quite liberal in its rejection of redundant search-tree branches. A detailed description of these subprocedures can be found in Appendix B.

However, we single out Append-Non-Redundant-Instances for discussion here, as this subprocedure contains the key logic which prevents factorial growth of intermediate results. We present this subprocedure in Algorithm 4.2. The rejection functionality is achieved with a simple check in lines 5 through 11 which ensures that within the current Least-Value-Set, only one representative from each symmetry subset is kept. This is done via an array Visited-Subsets which marks whether a subset has been visited; in lines 6 and 7, we use the absolute value of $\textsf{{{Symmetric-Subsets}}}\llbracket p\rrbracket$ as the index into Visited-Subsets. Note that for any $(p,q)$ within a Least-Value-Set, all of entries $\textsf{{{Values}}}\llbracket q\rrbracket$ in the Values array must be the same (since they are the least value), and thus any two entries belonging to the same symmetric subset must in fact be exchangeable (even if there may be labels with other values in the other slots of this symmetric subset!). In order to rely on this complete exchange equivalence (i.e., the slots are totally (anti)symmetric, and we likewise treat the labels as totally (anti)symmetric), we must use the symmetry-propagation mechanism of Section 4.1.1. And since we will now throw out any other intermediate configurations resulting from these symmetries, we will not be able to detect a future sign conflict, so it must have already been dealt with (via Zero-Due-to-Propagated-Symmetries).

4.4 Comments on complexity

Rather than perform detailed complexity analysis, we will present data in Section 5 which shows that our algorithm, like the original Butler-Portugal algorithm, is polynomial in the most common situations. Here we will make only a few comments which relate to our improvements to the algorithm.We have made two main improvements: First, we have chosen a means to represent the label symmetry group $L$ by distributing this information among several small arrays: Values, Label-Groups, and Propagated-Symmetries for storing the additional label symmetries which might be added during the course of the algorithm. Each of these arrays is just a single row of length $n$, and thus we incur an $O(n)$ memory cost, which is dwarfed by the $O(n^{3})$ cost of storing the Schreier-tree structure of the slot symmetry group $S$. While the new data structures require special logic to deal with them, their benefit is to reduce all group-theory computations involving the label symmetry group to either to array lookups which take $O(1)$ time, or single scans through the array which take $O(n)$ time. In either case, the time taken to compute with the label symmetry group is now insignificant.Second, the two arrays Symmetric-Subsets and Propagated-Symmetries allow us to eliminate equivalent choices from the search tree when presented with several interchangeable labels in a set of slots which are (anti)symmetric. At a cost of $O(n)$ space, we can avoid having to store $O(n!)$ intermediate results, which take $O(n!)$ time to process, a dramatic improvement.However, we note that $O(n!)$ behavior has not been eliminated entirely. Our algorithm is only capable of detecting subsets of indices which are totally symmetric or antisymmetric. But there are other possibilities which lead to combinatorial explosion. For example, suppose that a tensor has many indices which are pairwise symmetric, with a slot symmetry group generated by the pairwise exchanges

[TABLE]

Such a slot symmetry group has order $n!$ and will not be detected because it does not involve the direct exchange between two slots. One can imagine other factorial-size groups such as alternating groups or groups which symmetrize length-$k$ subsets, etc. These are not caught by the algorithm, and thus it will exhibit $O(n!)$ behavior in these cases. Fortunately, however, these cases are rare in actual calculation—in contrast to the (anti)symmetric case which we have addressed!One can see, in fact, that it is impossible to eliminate all possible sources of $O(n!)$ behavior (as should be expected [9]). Take the simple case of symmetric exchanges of length-$k$ subsets, for example. It may be possible, with some cleverly-designed data structures, to write an algorithm which—given data about these exchange symmetries ahead of time—can resolve canonicalization problems in polynomial time for such groups. But even if such an algorithm were designed, this is not enough, because that data about the symmetric exchange of length-$k$ subsets must first be somehow obtained. To generate it algorithmically for a specific $k$ requires $O(n^{k+1})$ group membership tests (and thus, in our particular case which resolves the question for length-1 exchanges, one can construct the Symmetric-Subsets structure using $O(n^{2})$ group-membership tests). To do so for all $k$ (thus, up to $k=n/2$) would then require $O(n!)$ group membership tests, and so we gain nothing; removing one head from the hydra only springs forth more.We point out, however, that there is a case in which data about length-$k$ exchanges can be provided ahead of time, in polynomial time: The case where a tensor monomial is constructed out of the tensor product of several identical factors,

[TABLE]

In such cases, it may be useful to have an algorithm which efficiently handles this type of symmetry. But even so, it will still be a special case, and $O(n!)$ behavior is still possible by some other means.

5 Performance testing

To test the effectiveness of Algorithm 4.1, we have fully implemented it and tested it in several situations against two implementations of the Butler-Portugal algorithm: the one in xPerm [4] which is used to simplify tensor expressions in xAct [5], and the built-in function TensorReduce which was introduced in Mathematica version 9. For the purpose of these tests, Algorithm 4.1 has been implemented in Mathematica code, and wrapped within the command Compile[…, CompilationTarget$\phantom{}\to\phantom{}$"C"]. This causes Mathematica to automatically generate C code (from the enclosed Mathematica code), which is then compiled to machine language via an external C compiler.222222It is therefore conceivable that one might achieve slightly better performance by coding directly in C, which would allow more control over memory resources. We have used the most up-to-date version of xPerm, from xAct version 1.1.2. The function TensorReduce is that available in Mathematica version 11.0.0.0, which was used to run all of these tests.In each test, we run each of the three implementations on randomized input with certain constraints as will be explained below. We have attempted to increase the input size until the asymptotic behavior becomes clearly visible, or at least until the total running time begins to exceed a few minutes. For each input size, we run the test several times (typically 10, except for the Riemann tensor tests which used 30 trials in order to generate enough results for each of the zero and nonzero bins). Different random input is generated for each trial. In each graph, the error bars represent the full range of running times over the trials for a given input size. The plot marks are placed at the geometric mean of the running times (thus, in these log plots, representing the mean of the logs).Each of the three implementations takes slightly different input, but we have made our best effort to eliminate as much overhead as possible. Algorithm 4.1 takes a pre-computed Scheier tree structure as input, as well as Symmetric-Subsets, and represents the label symmetry group via the inputs Values and Label-Groups. For xPerm, we have given it a pre-computed base-and-strong-generating-set (BSGS) for the slot symmetry group, and checked a flag indicating that this does not need to be re-computed. Unfortunately with Mathematica’s built-in function TensorReduce, there is no way to specify that the slot symmetries are already a BSGS, and so it spends time running the Schreier-Sims algorithm attempting to construct a BSGS which dominates the timing for low numbers of slots. Therefore we should point out that the data here are not entirely “fair”; however, the differences in performance are quite significant, in some cases differing by a whole power of $n$ and not just a constant factor.We also point out that each of these three implementations produces slightly different output; that is, each algorithm is using a slightly different notion of what “canonical” means. xPerm works in two stages, first moving all of the free indices as far forward as possible, then fixing them in place and canonicalizing the remaining indices by lexicographic order. Mathematica uses a definition of “canonical” based on the slot numbers of the contracted pairs, rather than the labels that appear in those slots, and so it may produce a different ordering. Algorithm 4.1 uses similar criteria for “canonical” as does xPerm, but works in one stage only, moving free and dummy indices as they appear; it is possible for them to give different results. Nevertheless, we have tested Algorithm 4.1 extensively for correctness, before running these timing tests.The Mathematica source code used to implement the algorithm and run these tests can be found at https://github.com/bniehoff/tensor-canonicalizer.

5.1 Preliminaries

We start with three tests which do not significantly stress any of the algorithms, and give a sort of basis of comparison which may help counteract the imbalance just mentioned. First, in Figure 5, we run the algorithms with free indices, with total slot symmetry. Thus there is no label renaming ambiguity; the algorithms need only sort the list of indices. From Figure 5 it is clear that Mathematica is doing extra work by running Schreier-Sims every time; furthermore it is clear that at around 50 indices, Mathematica changes strategy to use a randomized Schreier-Sims algorithm. To obtain an estimate on the asymptotic complexity, we attempt to fit a power law to the data points with at least 100 indices (thus ignoring the transient effects for low numbers of indices). We obtain the following power law growth:

[TABLE]

The slightly smaller power in the growth rate for Algorithm 4.1 compared to Mathematica is likely due to the fact that Algorithm 4.1 starts from slot-symmetry data where the orbits $\Delta^{S}_{i}$ are pre-computed (we point out that in a system where tensors are declared in advance, it is always possible to do this). Another possible explanation is that the original version of Butler-Portugal canonicalizes the free indices separately from the dummy ones using a slightly different algorithm [4, 30].

Next we run two basic tests with dummy indices in Figure 6. In the first test there are no slot symmetries, but instead we give a random sequence of dummy indices. Thus, all the dummies need to be renamed in the optimal order, but since there are no slot symmetries, there should be no ambiguity in the search tree. In Figure 6(a), one can clearly see the advantage of our data structures Values and Label-Groups for representing the label group $L$, which allow us to reduce the computation of the label-group orbit $\Delta^{L}_{g\llbracket j\rrbracket}$ to a simple array lookup. Using a similar method of best-fit power law while throwing out the points with a small number of indices, we obtain the time complexities

[TABLE]

In Figure 6(b) we run a test which eliminates both the advantages of Algorithm 4.1 and the drawbacks of the original Butler-Portugal algorithm. The algorithms are given the task to canonicalize the expression

[TABLE]

where the tensors $T,U$ have a cyclic symmetry over all of their indices. In the Mathematica algorithm, since the strong-generating-set of a cyclic group only has one generator, the Schreier-Sims step finishes quickly without significantly impacting performance. In Algorithm 4.1, one gets no benefit from the symmetry-propagation mechanism, because there are no (anti)symmetric subsets of slots. Nevertheless, the cyclic symmetry, combined with dummy renaming, provide enough computational challenge that it dominates in both cases, and Figure 6(b) shows nearly identical performance.232323The behavior of xPerm in Figure 6(b) is somewhat strange in that there is a sharp corner around 35 dummy pairs where it suddenly jumps to a different power-law behavior. We have discussed this with the author of [4] but cannot ascertain the cause. xPerm does not show any strange behavior in any other tests. The best-fit power law complexity is around $O(n^{3})$.

5.2 Random Riemann invariants

The next test is to repeat the benchmarking done in [10, 4, 31] for randomly-contracted Riemann monomials of the form

[TABLE]

where each $\circ$ should be filled with a random dummy index.242424Unlike [31], we do not include derivatives of the Riemann tensor, but only algebraic invariants. We also do not impose any dimension-dependent identities, but use only the Riemann tensor symmetries. Since the Riemann tensor has antisymmetries, these invariants are occasionally zero. The cases which turn out to be zero tend to have vastly shorter timing, and so we separate results into two bins, as shown in Figure 7. In Figure 7(a), for nonzero results we obtain the best-fit complexities

[TABLE]

We observe that not only is Algorithm 4.1 faster than Mathematica by at least an order of magnitude, it is also faster by about $n^{1.5}$ in time complexity. It appears xPerm is not well-fit by a power law in this case, although perhaps we have only managed to capture its transient behavior.For invariants which turn out to be zero, Figure 7(b) shows an even more dramatic improvement of up to 3 orders of magnitude and the following best-fit complexities:

[TABLE]

All of the algorithms are faster; again we see that Algorithm 4.1 is faster than Mathematica by about one power of $n$ and a large constant factor. Here we see the benefits of the symmetry-propagation mechanism for early zero detection.

5.3 Total symmetry

Next we run a test which exhibits the condition Algorithm 4.1 was specifically designed for: namely, contractions between totally symmetric tensors $T$ and $U$:

[TABLE]

as shown in Figure 8. We run two different versions of the test. First in Figure 8(a) is the “maximally frustrated” situation, where the dummy pairs are distributed so that each pair has one leg on $T$ and one leg on $U$:

[TABLE]

This situation represents the absolute worst-case scenario for the Butler-Portugal algorithm. In Figure 8(b), we do instead a fully randomized test, where the indices are distributed across both tensor factors randomly. This is a slightly better situation for Mathematica and xPerm, but the results are not much different. As expected, Algorithm 4.1 clearly exhibits better performance in this test, showing polynomial growth whereas Mathematica and xPerm grow as $O(n!)$. The data points for Algorithm 4.1 in the frustrated case of Figure 8(a) are well-fit by the power law

[TABLE]

While Mathematica and xPerm are clearly factorial, there are too few data points to get a good fit.

5.4 Pairwise symmetry

In the final test, we try a situation which is not addressed by our improvements, and in which we expect factorial growth in all cases. This test uses tensors $T$ and $U$ which have total symmetry among pairs of indices,

[TABLE]

and similarly for $U$. The results of this test are shown in Figure 9, again separated into the “frustrated” case in Figure 9(a) where each dummy pair is split between $T$ and $U$, and the randomized case in Figure 9(b) where the dummy pairs are randomly distributed across both tensors.

As expected, Figure 9(a) shows clear factorial growth, and is also the only case in these tests where Algorithm 4.1 performs worse than Mathematica, clearly growing at a faster rate. We expect that this is due to the extra processing done with Propagated-Symmetries which is wasted effort in this particular test case. By contrast, the performance for all three algorithms in Figure 9(b) is more equal (because when the labels are randomly distributed, not every pair of slots is necessarily equivalent to every other pair); nevertheless, for high numbers of indices we begin to see the patterns of Figure 9(a) repeated.We suspect that some adjustments in implementation could bring the performance of Algorithm 4.1 more in line with Mathematica in the pairwise-symmetric case (and similar) without sacrificing the performance benefits of the prior tests. However, we point out that real-world occurrences of the pairwise symmetry tested in Figure 9 are somewhat rare. If one did expect to see them frequently in some given application, it would be worth developing some more sophisticated version of the symmetry-propagation mechanism to bring them into the realm of polynomial time.

6 Discussion

It is important in many applications of symbolic computation to have algorithms which can efficiently deal with expressions with tensor indices. The Butler-Portugal algorithm [10] has been the state-of-the-art for some time, finding implementations in several popular software packages [4, 5, 7, 6, 3, 2]. We have presented two improvements on the Butler-Portugal algorithm: First, a faster way to deal with the label-exchange symmetries implicitly rather than explicitly; and second, a way to canonicalize contractions over subsets of symmetric or antisymmetric slots in polynomial time (where previously, they took factorial time).The improvements are evident in the data of Section 5 for all of the most common situations actually found in tensor calculations. These improvements are achieved using a number of different data structures: First we introduce the Values and Label-Groups arrays which not only allow an efficient specification of the problem (i.e., the partitioning of labels into free indices, component indices, and various separate pools of dummies of different types), but also allow one to reduce most label-group computations to simple array lookups. We also introduce the Symmetric-Subsets and Propagated-Symmetries arrays which allow efficient tracking of the subsets of slots which are totally (anti)symmetric, and how those symmetries interact with other slots in the expression via propagation over contracted dummies.In particular, the contractions-over-symmetric-subsets problem is dramatically reduced from factorial to polynomial complexity in both time and memory. This is achieved via the concept of symmetry propagation, wherein the slot symmetries at one end of a contraction can “travel” to the other end, where they become label symmetries. In this work, we have focused on propagating simple (anti)symmetries (i.e., those which are generated by the simple exchange of two slots), because these symmetries are the simplest to describe, and they address the worst-case situations likely to be encountered in real computations.However, we point out that the symmetry-propagation concept applies to generic symmetries as well. We strongly suggest this approach as an avenue to future improvement. The main challenge is designing an efficient means to describe such generic symmetries. One needs a data structure which is both simple to construct and simple to work with, and thus there are trade-offs to consider. We expect, however, that such an approach will be useful for cases like the pairwise symmetry of Section 5.4, provided that the pairs can be provided ahead of time (which is true, if these symmetries derive from the tensor product of many identical factors). We also expect that one may find efficient algorithms to enumerate Young-tableau-like symmetries, where one has totally (anti)symmetric exchange of totally (anti)symmetric length-$k$ subsets. Despite these potential future improvements, it is still true that the generic worst-case is $O(n!)$. However, if it is known in advance how such behavior could arise in a given real-world calculation, we are confident that the ideas presented here will be useful in developing further modifications that can reduce such behavior to polynomial.One should also consider the broader context in which tensor index canonicalization occurs: namely, to simplify algebraic expressions with tensor contractions. As we have pointed out, there are other approaches to this problem, based on matching terms rather then canonicalizing them [8], which may be faster in certain cases. Also, there is the largely-unexplored arena of multi-term symmetries, for which a graph-theory approach seems to be useful [12]. A complete tensor computer algebra system should probably combine aspects of all of these approaches.We look forward to addressing many of these considerations in future work.

Acknowledgements

I would like to thank J. M. Martín García for useful discussions and extensive comments on the initial draft. This work is supported in part by ERC grant ERC-2013-CoG 616732-HoloQosmos.

Appendices

Appendix A Efficient computation with permutation groups

A permutation group on a sequence of $n$ symbols can in principle be quite large, and one needs an efficient method to manage computations. In what has now become the standard method [9, 13, 14], one uses the notion of a strong generating set and its corresponding base. This is essentially a way to give a hierarchical structure to a group, which paves the way for doing many (though not all) problems of finite group theory in polynomial time.To define a base-and-strong-generating-set (or BSGS), let us first define a few group-theory notions. Suppose $G\subseteq S_{n}$ is a permutation group acting on a set $\Omega$ of $n$ symbols. Let $\alpha=\langle{\alpha_{1},\alpha_{2},\ldots,\alpha_{n}}\rangle$ be a fixed ordering of the points of $\Omega$. Then for example, a permutation $g\in G$ acts by rearranging the ordering to $\alpha^{g}=\langle{\alpha_{1}^{g},\alpha_{2}^{g},\ldots,\alpha_{n}^{g}}\rangle$. If any of the $\alpha_{i}$ are not moved by $g$, then we say $g$ stabilizes these $\alpha_{i}$. Given a subset of points $\{{\alpha_{i}}\}\subseteq\Omega$, the set of all $g\in G$ which stabilize these points give a subgroup of $G$ called the pointwise stabilizer of $\{{\alpha_{i}}\}$. The standard way to notate a pointwise stabilizer is with a subscript,

[TABLE]

which means that $G_{\alpha_{1}\alpha_{2}\ldots\alpha_{k}}$ is the pointwise stabilizer of all the $\alpha_{i}$ listed.Given the ordering $\alpha$, we can then consider the pointwise stabilizers of the subsequences of $\alpha$: $G_{\alpha_{1}}$, $G_{\alpha_{1}\alpha_{2}}$, $G_{\alpha_{1}\alpha_{2}\alpha_{3}}$, etc. Introducing a special notation for the stabilizer of the first $i-1$ points of $\alpha$,

[TABLE]

we can then observe that each $G^{(i)}$ contains $G^{(i+1)}$ as a subgroup, and thus the group $G$ can be decomposed into the point stabilizer chain

[TABLE]

The stabilizer chain gives a hierarchical way to organize a group which permits many efficient methods of computation.Given an ordering $\alpha$ and its corresponding stabilizer chain, the base of the group $G$ with respect to the ordering $\alpha$ is defined as the shortest subsequence of $\alpha$ whose pointwise stabilizer is trivial. That is, if $k\leq n$ is the lowest index for which $G^{(k)}=\{{\mathrm{id}}\}$, then the base with respect to this ordering is $B=\langle{\alpha_{1},\alpha_{2},\ldots,\alpha_{k-1}}\rangle$. Since the pointwise stabilizer of $B$ is the identity, each member $g$ of the group $G$ is uniquely identified by the image $B^{g}=\langle{\alpha_{1}^{g},\alpha_{2}^{g},\ldots,\alpha_{k-1}^{g}}\rangle$. Thus a base serves is a sort of index or set of coordinates.A strong generating set of $G$ with respect to a base $B$ is a set of permutations $S=\{{g_{1},g_{2},\ldots,g_{N}}\}$ which generate the group $G$, with an added condition that certain subsets of $S$ also generate each of the subsequence stabilizers $G^{(i)}$. Specifically, $S$ satisfies the property252525Here the notation $\langle\cdot\rangle$ means the group generated by the enclosed set. We trust this does not cause confusion with the usual meaning of $\langle{\ldots}\rangle$ in this article.

[TABLE]

It is important to note that, as a generating set for $G$, $S$ may have redundant generators which are needed in order to satisfy the property (A.4). For example, the symmetric group $S_{5}$ may be generated by the set (in cyclic notation)

[TABLE]

but this is not a strong generating set. To obtain the property (A.4) for the base $B=[1,2,3,4,5]$, we can use instead the set

[TABLE]

When constructing a strong generating set, it is still desirable to eliminate redundant generators, and there are standard algorithms for doing this.Given a base and a strong generating set, there are many computations that can be done efficiently. For example, one can test for group membership in $G$, find the order of $G$, or generate a random element of $G$ from a uniform distribution. To construct a BSGS in the first place requires the Schreier-Sims algorithm, whose main variant takes $O(n^{5})$ time. For detailed descriptions of these and other group theory algorithms, we refer the reader to [9, 13, 14].

Appendix B Subprocedures of the improved algorithm

Here we collect the various subprocedures of Algorithm 4.1 whose explicit code does not yield any great insights worthy of discussion in the main text, yet nevertheless whose function is not entirely obvious and ought to be written out.The first of these is Algorithm B.1 for the subprocedure Get-Least-Value-Instances, which traverses the current level of the search tree and generates the next level, by populating the Config-Slots-of-Least-Value array with all of the configurations $(g,s)$ and slot pairs $(p,q)$ at which one can find the least value label. In particular, in lines 8 through 16, we see how the slot number $q$ is obtained, by starting from $p$ and traversing any even-numbered symmetric sets available in Propagated-Symmetries. We only use the even-numbered sets because those ones are the result of symmetry propagation; the odd-numbered sets are already available in the slot symmetry group $S$, and thus should have been reached as one of the slots $p\in\Delta^{S}_{i}$.Next, in Algorithm B.2, we give the subprocedure Update-Propagated-Symmetries, which updates the Propagated-Symmetries array to reflect any new symmetric subsets found thus far in the search tree. Some effort must be taken to maintain the conditions of Definition 4.2.4 which ensure the correct behavior of Algorithm 4.2. The function Get-Next-Odd-Number in line 10 must maintain an internal state, so that it returns a new odd number every time it is called, in sequence, starting from 1. The function Remove-Singletons in line 27 removes any nonzero entry from Next-Propagated-Symmetries that occurs only once, thus preserving Rule P6 of Definition 4.2.4. Alternatively, one could avoid entering singletons in the first place, although this requires a bit more logic.In Algorithm B.3 we give the subprocedure Zero-Due-to-Propagated-Symmetries which looks for conflicts between Symmetric-Subsets and Propagated-Symmetries which may cause the expression to be zero. In line 12 is a check whether a propagated slot symmetry conflicts with an original slot symmetry. In line 21 is a check for dummy contractions where both ends lie in the same symmetric subset of slots, to see whether the metric (anti)symmetry (if there is a metric) agrees with the slot (anti)symmetry. We note that the functionality of Zero-Due-to-Propagated-Symmetries could also have been integrated into Update-Propagated-Symmetries, so that zero checks are made as soon as new symmetries are entered; this is the method used in the actual implementation tested in Section 5, although the logic becomes more complex.Finally, in Algorithm B.4 we give Update-Values-and-Label-Groups, whose purpose is to maintain the conditions of Definition 4.2.1 and Table 1 within the Values and Label-Groups arrays. Whenever the least-value label is consumed, we mark this by increasing the value of every label after it, within the set of labels whose value was least-value. This creates a new subset of labels with value $(\textsf{{{least-value}}}+1)$. The details of this behavior must be modified depending on whether the label is a component index, or a dummy index with or without metric. For dummy indices, we must offset the two entries at least-value and $(\textsf{{{least-value}}}+1)$, increasing those remaining to $(\textsf{{{least-value}}}+2)$. For example, take the Values array of (4.13):

[TABLE]

If we place the label $c_{1}$ into a slot, then we must update the Values array to

[TABLE]

This effectively marks the partner of a previously-used dummy for preferential treatment (as it must automatically precede the rest of the dummies in lexicographical order). In this case $c_{2}$ has a lower value than the remaining dummies $\{{d_{1},d_{2},e_{1},e_{2},f_{1},f_{2}}\}$, because $c_{2}$ is no longer interchangeable with them, and it should be preferred as it comes earlier in order.