Tensor network complexity of multilinear maps
Per Austrin Petteri Kaski Kaie Kubjas
cs.CC cs.DS





plainurl


































[name=Theorem,numberwithin=section]Thm
[name=Corollary,numberlike=Thm]Cor
[name=Lemma,numberlike=Thm]Lem
[name=Proposition,numberlike=Thm]Prop
[name=Claim,numberlike=Thm]Claim
[name=Conjecture,numberlike=Thm]Conjecture


Proof










Tensor network complexity of multilinear maps



School of Computer Science and Communication, KTH Royal Institute of Technologyaustrin@kth.se

Department of Computer Science, Aalto Universitypetteri.kaski@aalto.fi

Department of Mathematics and Systems Analysis, Aalto University
Laboratoire d'Informatique de Paris 6, Sorbonne Universitékaie.kubjas@aalto.fi

P. Austrin, P. Kaski, and K. Kubjas



Per Austrin, Petteri Kaski, and Kaie Kubjas



Theory of computation → Models of computation;
Theory of computation → Computational complexity and cryptography;
Theory of computation → Design and analysis of algorithms





An extended abstract of this paper appears in Proceedings of ITCS 2019.




Per Austrin was funded by the Swedish Research Council, Grant 621-2012-4546.  Petteri Kaski has received funding from the European Research Council, Grant 338077. Kaie Kubjas was supported by Marie Skłodowska-Curie  Grant 748354.


We are grateful to Andreas Björklund for highlighting branchwidth to us as a natural parameter to generalize from clique-counting to counting homomorphisms.







Tensor network complexity of multilinear maps
    Kaie Kubjas
    December 30, 2023
=============================================






  We study tensor networks as a model of arithmetic computation for
  evaluating multilinear maps.
  These capture any algorithm based on low border rank tensor decompositions, such as O(n^ω+ϵ) time matrix multiplication, and in addition many other algorithms such as O(n log n) time discrete Fourier transform and O^*(2^n) time for computing the permanent of a matrix.
  However tensor networks sometimes yield faster algorithms than those
  that follow from low-rank decompositions.  For instance the
  fastest known O(n^(ω +ϵ)t) time algorithms for
  counting 3t-cliques can be implemented with tensor networks, even
  though the underlying tensor has border rank n^3t for all t ≥ 2.
  For counting homomorphisms of a general pattern graph P into a host graph on n vertices we obtain an upper bound of O(n^(ω+ϵ)(P)/2) where (P) is the branchwidth of P. This essentially matches the bound for counting cliques, and yields small improvements over previous algorithms for many choices of P.

  While powerful, the model still has limitations, and we are able to
  show a number of unconditional lower bounds for various multilinear maps,
  including:
  
  
  (a) an Ω(n^(P)) time lower bound for counting homomorphisms from P to an n-vertex graph, matching the upper bound if ω = 2.
    In particular for P a v-clique this yields an Ω(n^⌈ 2v/3 ⌉) time lower bound for counting v-cliques, and for P a k-uniform v-hyperclique we obtain an Ω(n^v) time lower bound for k ≥ 3, ruling out tensor networks as an approach to obtaining non-trivial algorithms for hyperclique counting and the Max-3-CSP problem.

  
  (b) an Ω(2^0.918n) time lower bound for the permanent
    of an n × n matrix.

  






§ INTRODUCTION


One of the cornerstones of the theory of computation is the study of efficient
algorithms:

For a function f, how much time is required to evaluate f on given inputs?

Answering this question for almost any specific f is well beyond
reach of contemporary tools. For example, it is theoretically possible
that canonical NP-complete problems,
such as the Circuit Satisfiability problem, can be solved in linear
time whereas they are widely believed to require super-polynomial
(or somewhat less widely, exponential) time <cit.>.
The main reason for this barrier to quantitative understanding is that it is
very hard to prove lower bounds for explicit functions in general models
of computation such as circuits or Turing machines.

This situation withstanding, a more modest program to advance our
understanding of computation is to study restricted models of computation
that for many f of interest are simultaneously


  * general enough to capture the fastest-known algorithms for
f, and

  * restricted enough to admit proofs of strong unconditional time
lower bounds for f.

There is a substantial body of existing work that fits under this program,
ranging from the study of low-depth or otherwise restricted circuits (see e.g. <cit.>, Ch. 14)
to models of algorithm-design principles such as greedy algorithms,
backtracking, or dynamic programming
<cit.>, to linear or semidefinite
programming relaxations for hard combinatorial optimization problems
<cit.>.



 §.§ Multilinear maps


One class of functions f that are of substantial importance is the family of ℓ-linear maps (multilinear maps) from ℓ input vector spaces to an output vector space.
[Multilinear maps with ℓ=1 are called linear, ℓ=2 bilinear, ℓ=3 trilinear, and so forth.]
Examples range from maps of known near-linear-time complexity in the input size, such as the discrete Fourier transform <cit.>, to maps without known polynomial-time-complexity algorithms, such as the permanent of a matrix <cit.>. Beyond motivation in numerical multilinear algebra and its applications, recent advances in the study of fine-grained algorithm design and complexity have highlighted the fundamental role of algebraic methods in the fastest-known algorithm designs across a broad range of tasks from graph problems, such as all-pairs shortest paths and k-clique, to parsing and constraint satisfaction problems, such as maximum satisfiability and graph coloring <cit.>.

In this paper, we study the arithmetic complexity of evaluating
a multilinear map, that is, the number of operations (scalar additions,
subtractions, and multiplications) needed to evaluate the map.
To set up a baseline, a generic ℓ-linear map from ℓ vector spaces
of dimension n to a scalar requires Ω(n^ℓ) scalars to represent
the map directly using combinations of basis vectors. Given this
complexity of a direct explicit representation, it is a fundamental
problem to seek less costly representations, along with associated
efficient algorithms that work on the chosen representation.

We propose the systematic study of tensor networks on hypergraphs as
a framework for fast evaluation of multilinear maps, and show a number of
upper and lower bounds on the computational power of tensor networks
in the spirit of (i) and (ii) above.




 §.§ Tensor networks




Tensor networks have a long and rich history which can be traced as
far back as 19-century studies in invariant theory due
to Cayley <cit.>, Clebsch <cit.>,
Clifford <cit.>, Sylvester <cit.>,
and Kempe <cit.>.
Tensor networks are extensively deployed in applications from pure mathematics
and theoretical physics <cit.> to computational physics and chemistry <cit.>.
In theoretical computer science,
they appear in various guises including, for example,
in the Holant framework <cit.>,
in the study of probabilistic graphical models <cit.>,
in the study of parallel programming <cit.>,
in the study of quantum computing <cit.>,
and in the study of arithmetic
complexity <cit.>.
Tensor contraction is also emerging as a basic
computational primitive in computer
hardware <cit.>.
(See <ref> for a detailed discussion.)
As the precise definitions are somewhat technical, let us start with
a few simple motivating examples and then state our results,
with the understanding that precise definitions appear
in <ref>.

In our setting, a tensor network is a hypergraph in which the vertices are
tensors and the hyperedges are called modes.
Each mode that is incident to a tensor defines a “dimension”
for indexing the entries of the tensor—for example, a matrix is a tensor
that is incident to two modes, one mode for the rows of the matrix, and
the other mode for the columns of the matrix. A network may be simplified
by a sequence of contractions, where each contraction takes a subset
of tensors and replaces it with a single tensor whose entries are obtained
as generalized inner products of the entries of the tensors being contracted.

As a concrete first example of these concepts, let us consider the task of
multiplying two matrices, A and B. More specifically, let
A be a matrix with rows indexed by mode i and
columns indexed by mode k, and let B be a matrix
with rows indexed by mode k and columns indexed by mode j.
We may represent the multiplication task as
the tensor network depicted on the left in (<ref>).
The result of contracting A and B is a new
matrix with rows indexed by i and columns indexed by j, where the
entry at each position (i,j) is ∑_k A_ik B_kj.  If the three
index sets all have size n, then computing A · B by contracting
them in such a direct manner uses Θ(n^3) operations.  To obtain
faster matrix multiplication, we can rewrite the bare-bones network
on the left in (<ref>) using a low-rank
decomposition of the matrix multiplication tensor.
For example, Strassen's decomposition <cit.>
of 2 × 2 matrix multiplication
can be represented using the second network in (<ref>).
Note that the index i used by A and the result has been separated into
two distinct indices i and i', and similarly for j and k.
0.9figures/intro_matrixmul1intro matrix mul 1

We can execute the network by succesively contracting
groups of vertices.  In (<ref>) we see the
process of successively contracting pairs of tensors in a carefully
chosen order, until only a single tensor – the result of the
computation – remains.  Such an execution can be naturally
represented by a rooted binary tree, as shown on the right in
(<ref>), where the tensors of the network form
the leaves, and each internal node represents the result of
contracting its two children.  To summarize, a tensor-network
algorithm is specified by providing (a) a tensor network that when
contracted yields the desired result, and (b) an execution tree
indicating the order in which to contract tensors in the network.

0.74figures/intro_matrixmul2intro matrix mul 2

The cost of performing one of the contractions in an execution
is the product of the lengths of the modes used by any tensor involved
in the contraction.  This simply measures (up to a constant factor)
the number of arithmetic operations (additions/multiplications) used
to compute the result by a direct, naïve computation that does not
depend on the values of the tensors.  For example, the contraction of
α and A in the first step of (<ref>)
has cost 28 because it involves the three modes i' (length 2),
k (length 2) and ℓ (length 7).

We observe that cost is data-oblivious – the tensor α is fixed with many
zero-entries but these entries still contribute to the cost.  Indeed, in many cases there may be faster ways of evaluating a contraction than to evaluate it naively, and just like we saw above, this can often be dealt with by rewriting the network appropriately.
For instance, consider now the multiplication of two 2^k × 2^k
matrices.  Because the family of matrix multiplication tensors is closed under
Kronecker products, this operation may be computed by a tensor
network like the one shown in (<ref>)
(depicting the case k=5), where α, β and γ are
as in (<ref>).  The rows/columns of the
matrices are now indexed by k-tuples of bits.  The execution of this
network contracts one α/β/γ tensor at a time,
which lets us keep the cost low.  For example, the first contraction of A
with the first α block has a cost of 2^k · 2^k · 7,
and results in a tensor of size 2^k-1× 2^k-1× 7,
then the next contraction has a cost of 2^k-1· 2^k-1·
7^2 and produces a result of size 2^k-2× 2^k-2× 7
× 7, and so on, until the contraction with the last α
block which has a cost of 2 · 2 · 7^k = O(7^k), and all the
contractions in the execution have cost bounded by this, meaning that
we get a total running time of O(k 7^k) = O(N^log_2 7log N) for
N × N matrices.
[In fact, a more careful analysis gives running time O(N^log_2 7).]
0.65figures/intro_matrixmul3intro matrix mul 3

This type of argument can capture any algorithm based on a low-rank
decomposition of the underlying tensor of the multilinear map, and
indeed, this enables O(n^ω)-time[Throughout the paper,
ω = ω() denotes the infimum over all e such that
the arithmetic complexity of multiplying two n × n matrices is
O(n^e).  While the value of ω may depend on the underlying
field , we tacitly ignore this, since the field is fixed throughout
the paper. For all fields it is known that
2≤ω<2.3728639 <cit.>.]
 matrix multiplication
using tensor networks.
Beyond simple low-rank decompositions, which always give rise to
“star-like” networks as in (<ref>), there are
many interesting algorithms that can be captured using networks with a
more complicated topology.
For instance, many maps of substantial
importance have a layered structure that decomposes the map to a sequence
of elementary maps. A canonical example is the
discrete Fourier transform (DFT),
which for a smooth composite order such as 2^k, can be decomposed
into a fast Fourier transform (FFT) that consists of a sequence
of k transforms of order 2 interleaved with diagonal-matrix
multiplications of twiddle factors <cit.>.



 §.§ Our results



Starting with motivation (i) and seeking to express existing fastest-known
algorithms
as executions of tensor networks by a sequence of contractions,
we show upper bounds for a number of
natural problems. Beyond standard linear settings such as the FFT,
not only do tensor networks realize classical bilinear settings
such as Abelian group algebra products and
fast matrix multiplication algorithms based on low tensor rank,
they are powerful enough to capture a substantial number of higher-linearity
applications, including Ryser's algorithm for matrix permanent <cit.>,
and the Kruskal operator <cit.>
(cf. <ref>),
which underlies realization of rank-decompositions for tensor rank
<cit.> and current fastest algorithms for detecting outlier
correlations <cit.>.

One problem for which tensor networks turn out to be particularly useful
is counting homomorphisms from a fixed pattern graph P to
a large host graph G on n vertices.  The most well-studied such problem
is when P is a k-clique.  For this problem, the currently fastest
algorithm runs in time roughly O(n^ω k/3) (with variations in
the exponent depending on k  3)
<cit.>.  For general P, it is known
that the problem can be solved in O(n^(P)+1) time
<cit.>, where (P) is the treewidth of P.
We show that tensor networks can solve the problem in O(n^(ω+ϵ)
  (P)/2) time, where (P) is the branchwidth of P.
For P a k-clique we have (P) = ⌈ 2k/3⌉ so this
almost recovers the O(n^ω k/3) running time, and in this case
we can slightly improve the running time to recover the
O(n^ω⌊ k/3 ⌋ + (k  3)) time of Nešetřil and
Poljak <cit.>.  In the case of general P, this improves
on the treewidth-based bound for graphs with (P) ≤
2((P)+1)/ω (and in particular if ω=2 it is always as
fast as the treewidth-based bound, ignoring the ϵ).  By
recent results of Curticapean, Dell, and Marx <cit.>,
fast algorithms for homomorphism-counting can be used to obtain fast
algorithms for counting subgraphs of G isomorphic to P, and in
some cases our new branchwidth-based bound leads to an improvement;
for example, for counting paths of lengths of length 7, 8 or 9,
we get a running time of O(n^3ω/2+ϵ) < O(n^3.56)
compared to O(n^4) using the treewidth-based bound, whereas for very long
paths it is not clear whether we would need ω=2 in order for
these bound to improve on the treewidth-based bound. Previous work
that combines branch decompositions and fast matrix
multiplication includes Dorn <cit.> and
Bodlaender et al. <cit.>.

Further applications
captured by tensor networks are the set covering and set partitioning
frameworks via fast zeta and Möbius transforms that underlie the
current fastest algorithms for graph coloring <cit.>
and its generalizations such as computing the
Tutte polynomial <cit.>.
To summarize, we have the following compendium theorem of upper bound results.
For the detailed definitions of the relevant multilinear maps,
cf. <ref> and <ref>.


  
  We have the following upper bounds on arithmetic complexity via tensor networks:
  
  
  * O(n^ω + ϵ) for the matrix multiplication map of two n × n matrices.
  
  * O(n^(ω+ϵ)⌊ v/3 ⌋  + (v  3)) for counting v-cliques in an n-vertex graph.
  
  * O(n^(ω+ϵ) (P)/2) for counting homomorphisms of a fixed pattern (hyper)graph P into a (hyper)graph on n vertices.
  
  * O(max(n^⌈ℓ/2⌉(ω+ϵ-1)r, n^2⌈ℓ/2⌉r^ω+ϵ-2)) for the Kruskal operator of ℓ matrices of shape n× r.
  
  * O(2^k k) for the discrete Fourier transforms for the Abelian groups _2^k and _2^k.
  
  * O(2^kk) for group algebra products on [_2^k] and [_2^k] when 2 is unit in .
  
  * O(2^kk) for the semigroup algebra product on [({0,1}^k,⊆,∩,∪)].
  
  * O(2^n n) for the permanent of an n × n matrix.

  
Above ϵ>0 is an arbitrary constant.


Perhaps the most interesting application above
is the v-clique problem, which suggests that one should seek to
pursue generalizations to v-vertex hypercliques of
vk hyperedges with v>k≥ 3.
Indeed, subgraph counting is a problem that has received substantial
interest over the years (e.g. <cit.>), but progress in the particular case of v-clique
has been stuck to the extent that the problem has attracted notoriety as
a hardness assumption in fine-grained
complexity <cit.>.
Beyond the study of cliques, hypercliques, and subgraph counting, nontrivial
algorithms for such forms would have immediate applicability, for example,
in the study of maximum constraint satisfaction problems (Max-CSP) for
constraints of width k≥ 3; cf. Williams <cit.> for
the case k=2. One of the main goals of our subsequent lower bounds is
to rule out tensor networks as a candidate to yield improved algorithms
in this setting.

Turning from upper bounds to lower bounds and motivation (ii), tensor
networks are restricted enough to enable nontrivial lower bounds for many
multilinear maps.
To begin with, an immediate limitation of tensor networks is
that all the intermediate results during the execution of a network are
multilinear, and the execution of a network can be simulated by a multilinear
circuit. Raz <cit.> shows that multilinear formulas cannot compute the
determinant of an n× n matrix in a polynomial number of operations
in n, even though polynomial-size general circuits are known for
the determinant (cf. <cit.>).

It turns out that considerably stronger lower bounds can be shown
for tensor networks. In particular, we establish essentially tight lower bounds
(subject to the assumption ω=2) for arithmetic complexity via tensor
networks of P-homomorphism counting and the Kruskal operator.
Furthermore, we rule out the possibility of any
nontrivial algorithm designs via tensor networks for counting cliques in hypergraphs.
The following theorem collects
our main lower-bound results, and should be contrasted with the upper bounds
in Theorem <ref>.


  
  We have the following lower bounds on arithmetic complexity via tensor networks:
  
  
  * Ω(n^(P)) for the multilinear form corresponding to P-homomorphism counting.  In particular, this yields a lower bound of Ω(n^⌈ 2v/3 ⌉) for counting cliques of size v, and a lower bound of Ω(n^v) for counting hypercliques of size v.
  
  * Ω(max(n^ℓ, n^⌈ l/2⌉r)) for the Kruskal operator of ℓ matrices of shape n× r.
  
  * Ω(nn/3) for the determinant or permanent of an n × n matrix.
  


We remark that <cit.> independently showed that the border rank of the v-hyperclique tensor is Ω(n^v); our Ω(n^v) lower bound for tensor networks strengthens that.
One may wonder about the gap between the bounds of
Ω(nn/3) and O(2^n n) for the permanent.
As we explain below, our lower bound methods are inherently rank-based
and cannot go beyond nn/3.  A curious point is that it is not immediately clear whether tensor networks can even achieve O^*(2^n) time for the determinant, and we do not know whether or not this is possible.




 §.§ Overview of proof ideas




As a running example in this overview, we consider the 6-linear forms
A: ^[n] × [n]×^[n] × [n]×…×^[n] × [n]→
taking as input 6 matrices of size n × n, defined by

    A(X^(1), X^(2), X^(3), X^(4), X^(5), X^(6)) = ∑_i, j, k, ℓ∈[n] X^(1)_ij X^(2)_ik X^(3)_iℓ X^(4)_jk X^(5)_jℓ X^(6)_kℓ .

If χ is the adjacency matrix of a loopless graph G, then A(χ, χ, χ, χ, χ,
χ) counts the number of 4-cliques in the graph.  Associated with A is the 6-tensor T(A) of size n^2 × n^2 ×⋯× n^2, where each of the 6 modes is indexed by a pair (i, j) ∈
[n] × [n], and the value at a given position is the coefficient of the corresponding term in A.  Concretely,

    T(A)_i_1j_1, i_2k_2, i_3ℓ_3, j_4k_4, j_5ℓ_5, k_6ℓ_6=
      1    if i_1 = i_2 = i_3, j_1 = j_4 = j_5, k_2 = k_4 = k_6∧ℓ_3 = ℓ_5 = ℓ_6,
    
      0    otherwise.




    Upper bounds


Most, but not all, of the families of multilinear maps we consider are
closed under taking Kronecker products.  For instance, consider the
4-clique counting form (<ref>) for an n-vertex
graph and its associated tensor T(A).  Then for any
k ≥ 1, the tensor associated with the 4-clique counting form in
n^k-vertex graphs is T(A)^⊗ k, the k-fold Kronecker
product of T(A) with itself.  We write A^⊗ k for the
associated map.  With this in mind, it is natural to seek general
constructions that, given an efficient evaluation
of some map A, yields an efficient evaluation of A^⊗ k.

We give such a construction, and show that the cost of the best tensor
network execution for A^⊗ k is essentially submultiplicative
in a quantity that we call the amortized cost of an execution.
For tensors of order at most 3, the notion of amortized cost
essentially captures the rank of T(A), but for higher-order tensors,
the amortized cost may be significantly smaller than the rank.
Roughly speaking, the amortized cost of a step in an execution of a
map A is: (i) equal to the normal cost if the operation involves the
contraction of two tensors that both depend on some input variables of
A, but (ii) equal to the size of the result if only one of the
tensors involved in the contraction depends on the input variables of
A.  A precise definition appears in <ref>.
Our general upper bound for the cost of A^⊗ k can, somewhat
informally, be stated as follows.

[Submultiplicativity of cost, informal statement]ThmGeneralUpperBoundTheorem
  
  If a multilinear map A has a tensor network execution
  consisting of s steps, each with cost at most c and amortized
  cost at most a, then A^⊗ k has a tensor network
  execution consisting of at most k · s steps, each with cost at
  most a^k-1· c.


An immediate corollary of this is that we can capture any algorithm
for A^⊗ k based on a low-rank decomposition of T(A)
(<ref>).  For example, this implies that
tensor networks can multiply n × n matrices in O(n^ω +
  ϵ) time (<ref>).

However, returning to our running example form (<ref>), as we explain below the tensor T(A) has rank n^4, meaning
that <ref> only yields a trivial upper bound.
This is where the full generality of <ref>
comes in.  Consider the form (<ref>) for graphs on
some constant number n_0 of vertices.  As it turns out, we can design
a network and an associated execution for this form, depicted
in (<ref>) and explained in more detail
in the proof of <ref>,
with an execution of cost n_0^2e+3 and amortized cost
n_0^e+1, where n_0^e is the rank of the tensor associated with
n_0 × n_0 matrix multiplication.  Picking n_0 to be a large
enough constant so that
e is approximately ω, and letting k be such that n is
approximately n_0^k, we obtain via <ref> an
O(n^ω +ϵ + 1) time upper bound for (<ref>).
0.75figures/intro_4choose2overview 4 choose 2 execution



    Lower bounds


Just like many other arithmetic complexity lower bounds, our lower bounds
boil down to establishing lower bounds on the rank of certain
matrices.

In order to establish a lower bound on the rank of T(A),
we flatten it to a matrix and analyze the rank of that
matrix.  There are 2^5 possible bipartitions of the 6 modes of
T(A) into two non-empty subsets, and the lower bound on the rank of
T(A) that we obtain is the maximum of the ranks of the
resulting matrices.  Using this method it is easy to establish that
for our example form (<ref>), the rank of
T(A) = n^4.  That this is an upper bound follows from
(<ref>), and that it is a lower bound follows by
considering the bipartition taking variables X^(1) and
X^(6) as row indices, and the other 4 variables as column indices.
The resulting n^4 × n^8 matrix has full rank.

Tensor networks are more versatile and can be more efficient
than low-rank decompositions of T(A).  Nevertheless, we show
limitations on this versatility.  In particular we show that every
tensor network execution for A induces a tree in which the leaves are the inputs of A and all internal vertices have degree 3.  We call this a socket tree.  Each edge in a socket
tree induces a bipartition of the variables and our key technical
lemma is to show that for each such bipartition, the rank of the
corresponding flattening of T(A) is a lower bound on the
cost of the execution that gave rise to the tree.  Thus, to obtain a lower
bound for the cost of a specific execution, we consider the maximum
rank obtained over all
edges of the corresponding socket tree, and to lower bound the cost of
every tensor network execution, we minimize this quantity over all
possible socket trees.  We refer to the resulting quantity as the
socket width of A, denoted w(A) (formal definition appears in <ref>).  Our general lower bound
can thus be phrased as follows, where c(A) denotes the minimum cost of a tensor network for evaluating A (formal definition appears in <ref>).

ThmGeneralLowerBoundTheorem
  
  For every multilinear map A, it holds that c(A) ≥ w(A).


Indeed, for our running example (<ref>), there are
low-width socket trees establishing that w(A) ≤ n^3, see
(<ref>).  However, that bound is tight: our
Ω(n^⌈ 2 · 4 / 3 ⌉) = Ω(n^3) lower bound
for the 42-linear form (<ref>)
is obtained by proving  that w(A) ≥ n^3 (<ref>) and
appealing to <ref>.
1.0figures/intro_4choose2_treeoverview socket tree







 §.§ Earlier and related work



We now proceed to a more detailed discussion of earlier work.



    Tensor networks


The history of tensor networks
(or, alternatively, tensor diagrams or diagrams) as
an analytical and computational tool goes back to the 19-century to the works
by Cayley <cit.>, Clebsch <cit.>,
Clifford <cit.>, Sylvester <cit.>,
and Kempe <cit.>. The diagrammatic form used here can be
traced back to Penrose <cit.>. Some early appearances of tensor
diagrams are by Cvitanovic̀ <cit.>,
Cvitanovic̀ and Kennedy <cit.>,
Kuperberg <cit.>, and many others. Surveys of tensor
diagrams can be found in Penrose and Rindler <cit.> and
Landsberg <cit.>. Schrijver <cit.> gives
a brief historical account.

A principal deviation from Penrose's notation that we make in this paper
is to work subject to a basis and a corresponding dual basis in each relevant
vector space to avoid distinction between primal and dual spaces.
This in particular enables a concise treatment of hyperedges and saves us
from considering orientation of edges, or the planar layout of edges
in a drawing. That is, we will view a tensor network combinatorially as a
hypergraph with tensors associated at the vertices, and with a subset
of the hyperedges designated to form the
boundary of the network (cf. <ref> for the
precise definitions). A yet further conceptual difference is that we view
the execution of a tensor network as a sequence of contractions
of sets of vertices (tensors) rather than as contractions of
hyperedges (modes). This choice enables us to reduce the size of hyperedges
gradually before eliminating them during an execution, thus
enabling better granularity. For simplicity, we will restrict to
a purely multilinear framework and will not consider sums of networks
although such a study would be possible, and is pursued e.g. in
Penrose's work <cit.>.

A large body of existing work in applications (cf. <ref>) studies
how to efficiently execute a given tensor network D.
Our quest in this paper differs from such studies in that
we study a multilinear map A, and seek to
design the most efficient network D that realizes A,
or to establish lower bounds for best-possible designs.
In particular, our upper and lower bounds in <ref>
and <ref> are over all tensor networks that realize
a particular map A of interest.



    Computational problems on tensors and tensor networks


Problems on given tensors and tensor networks are known to be mostly
computationally hard as soon as the setting changes from matrices to
higher-order tensors.
Håstad <cit.> showed that computing the rank of
an explicitly given ℓ-tensor is NP-complete over finite fields
and NP-hard over rationals whenever ℓ≥ 3. Hillar and Lim extended the latter result to any field containing ℚ <cit.>. They also showed that many other tensor problems such as the eigenvalue, singular value and spectral norm decision and approximation problems as well as rank-1 approximation problem for 3-tensors (over ℝ and in some cases over ℂ) are NP-hard.

The task of finding the value of a given scalar-valued tensor network
is known to be #P-complete (see e.g. <cit.>).
Similarly, it is NP-hard to find the most efficient sequence of contractions
for a given network <cit.>.



    Tensor networks in applications


Beyond our present use of tensor networks as a model of computation to efficiently evaluate multilinear maps, tensor networks are used across a broad range of applications. Accordingly, the following should be viewed as mere pointers to further literature on tensor networks, not as an exhaustive listing of all applications of tensor networks. Orus <cit.> gives an introduction to tensor networks in the context of computational physics. Itai and Landau <cit.> study quantum computation and quantum algorithms for evaluating tensor networks <cit.>. Solomonik and Hoefler study sparse tensor algebra as a model for parallel programming <cit.>. The Holant framework
introduced by Valiant <cit.> and studied further by
Cai et al. <cit.> involves the study of multilinear
sum–product expressions that can be viewed as tensor networks.
Tensor networks appear naturally in the study of probabilistic graphical
models <cit.>,
and in the study of various machine-learning
problems <cit.>.



    Bilinear and multilinear complexity


As was concisely outlined in <ref>, for bilinear maps our present work reduces to the study of tensor rank of 3-tensors and an extensive body of work on bilinear complexity, with the arithmetic complexity of the matrix multiplication map as the primary driver. For two starting points to this literature, we refer to the monograph by Bürgisser, Clausen, and Shokrollahi <cit.> and to the survey by Pan <cit.>. Our present work can be seen as generalizing this bilinear theory to higher orders of linearity via tensor networks and their executions. The current state of the art for fast matrix multiplication can be found in Le Gall <cit.>, Le Gall and Urrutia <cit.>, Vassilevska Williams <cit.>, Cohn and Umans <cit.>, Cohn, Kleinberg, Szegedy, and Umans <cit.>, Benson and Ballard <cit.>, and Huang, Rice, Matthews, and van de Geijn <cit.>.




 §.§ Organization of this paper


Section <ref> recalls preliminaries on tensors, multilinear maps, and branchwidth. Section <ref> reviews the specific multilinear maps that we study in this work and describes for each map its associated tensor. In Section <ref>, tensor networks, execution and cost of a tensor network, and cost of a multilinear map are defined. Section <ref> presents tensor-network algorithms for the multilinear maps introduced in Section <ref>. In Section <ref>, a general lower bound on the cost of evaluating a multilinear map using tensor network is obtained. The lower bound is expressed in terms of the socket-width of a multilinear map. In Section <ref>, lower bounds on socket-width for concrete multilinear maps studied in Sections <ref> and <ref> are derived. Finally, Appendix <ref> gives some background results on minimum-cost executions.



§ PRELIMINARIES



Throughout the paper [n] denotes {1,2,…,n} and  denotes
an arbitrary fixed field.



 §.§ Tensors



This section sets up our notation for tensors and multilinear maps.
We work with tensors and multilinear maps relative to fixed bases for
the respective vector spaces over .



    Modes, indexing, and positions


We will work with the following convention of positioning individual
entries inside a tensor.  Let E be a finite set of modes. Associate with each mode e∈ E a finite nonempty index set J(e). In this case we say that E is a set of indexed modes.  The length of e is |J(e)|. A position
  is an element j∈∏_e∈
  EJ(e). Let us write J(E)=∏_e∈ EJ(e) for the set of all
positions with respect to the indexed modes E.  In the special case
that the set of modes E is empty we define the set of positions
J(E) to consist of a single element.



    Tensors, matrices, vectors, and scalars




Let E be a set of indexed modes. A tensor with respect to E is
a mapping T:J(E)→. Equivalently, we write T∈^J(E)
to indicate that T is a tensor with respect to the indexed modes E.
We view the set ^J(E) of all tensors over E as a vector space over 
with addition and scalar multiplication of tensors defined entrywise.

We say that |E| is the order of the tensor.  A tensor of order
zero is called a scalar, a tensor of order one is called a vector, and a tensor of order two is called a matrix.

The volume of the tensor is |J(E)|. The tuple (|J(e)|:e∈ E)
is the shape of the tensor. It is convenient to use the
“×”-symbol to punctuate the shape of a tensor; that is,
instead of writing, say (2,3,4) for the shape, we write 2× 3× 4.
For a position j∈ J(E) and a tensor T∈^J(E),
we say that T_j∈ is the entry of T at j.

A flattening of T induced by a bipartition E_1 ∪ E_2 = E
of the modes of T is a |J(E_1)| × |J(E_2)| matrix M where, for j_1 ∈ J(E_1) and j_2 ∈ J(E_2) we have M_j_1,j_2 = T_j_1j_2.
Given two order ℓ tensors S ∈^[n_1] × [n_2] ×⋯× [n_ℓ] and T ∈^[m_1] × [m_2] ×⋯× [m_ℓ], their Kronecker product S ⊗ T is a tensor in ^[n_1 m_1] × [n_2 m_2] ×⋯× [n_ℓ m_ℓ] defined by

    (S ⊗ T)_m_1(i_1-1)+j_1,m_2(i_2-1)+j_2,…,m_ℓ(i_ℓ-1)+j_ℓ = S_i_1,i_2,…,i_ℓ T_j_1,j_2,…,j_ℓ.




 §.§ Multilinear maps


Let E_1,E_2,…,E_ℓ,E' be pairwise disjoint sets of indexed modes
such that E_1,E_2,…,E_ℓ are nonempty. We say that a map
A:^J(E_1)×^J(E_2)×⋯×^J(E_ℓ)→^J(E')
is an ℓ-linear map if
A is linear with respect to each of its ℓ inputs
individually when the other inputs remain fixed.
In particular, a 1-linear map is a linear map.
A multilinear map that gives a scalar output is a multilinear form.
In particular, A is a form if and only if E' is empty.



    The tensors of a multilinear map


For an ℓ-linear map
A:^J(E_1)×^J(E_2)×⋯×^J(E_ℓ)→^J(E') ,
we define two slightly different tensors T(A) and T̂(A).  Both are indexed by J(E_1 ∪ E_2 ∪…∪ E_ℓ∪ E') and at position j_1j_2… j_ℓ j' take the value

    T(A)_j_1j_2… j_ℓ j' = T̂(A)_j_1j_2 … j_ℓ j' =
    A(e^(j_1),e^(j_2),…,e^(j_ℓ))_j' ,

where e^(j_i)∈^J(E_i) denotes the tensor with a 1 in
position j_i and 0s in all other position.  The difference between
T(A) and T̂(A) is their shape.  The shape of T(A) is
|J(E_1)| × |J(E_2)| ×⋯× |J(E_ℓ)| ×
|J(E')|,
except if A is a form in which case the |J(E')| part is omitted.
Thus T(A) is of order ℓ+1 (or ℓ if A is a form).  The
shape of T̂(A) is (|J(e)|: e ∈ E_i, i ∈ [ℓ+1]), thus
its order is |E_1| + |E_2| + ⋯ + |E_ℓ| + |E'|.

In other words, each mode of T(A) corresponds to one of the inputs
of A, or the output.  These inputs are in turn sets of indexed modes
so may contain more “fine-grained” structure, but this information
is lost at the level of granularity of T(A).  When working with
tensor networks for evaluating A, we need to keep track of the
fine-grained mode structure because this is in many cases what allows
us to construct efficient algorithms, hence in most parts of the paper
we are more interested in the tensor T̂(A) which contains this
fine-grained structure.

On the other hand, T̂(A) does not contain information about
which modes are grouped together to form the inputs and output of A,
and this information is also important.  This leads us to the notion
of sockets, defined next.



    Sockets


Let us study the tensor T̂(A) with respect to the map A.
We say that the modes in
E_1∪ E_2∪⋯∪ E_ℓ are the input modes of T̂(A),
and the modes in E' are the output modes of T̂(A) with respect to A.
Let us say that E_1, …, E_ℓ are the
input sockets of T̂(A) with respect to A. Similarly,
E' is the output socket in T̂(A) with respect to A.
In particular, the output socket is empty if and only if A is a form. To describe a socketing of the modes of a tensor, it is convenient to
use parentheses to group a “×”-punctuated shape of a tensor
into sockets, see also <ref>.

Let T̂ be a tensor over a set of indexed modes E.
Any tuple ℰ = (E_1,E_2,…,E_ℓ,E') of subsets of E
that partitions E with E_1,E_2,…,E_ℓ nonempty
defines an ℓ-linear map
A_ℰ(T̂):^J(E_1)×^J(E_2)×⋯×^J(E_ℓ)→^J(E') with T̂(A_ℰ(T̂))=T̂. In this case the tuple
(E_1,E_2,…,E_ℓ) gives the input sockets of T and E'
is the output socket of T̂ with respect to A_ℰ(T̂).

We thus conclude that two multilinear maps A_1,A_2
may have the same base tensor T̂(A_1)=T̂(A_2), and from a tensor T̂ one may
obtain different multilinear maps by varying how the modes of T̂
are assigned to input and output sockets.



    The form of a multilinear map


Let A be a multilinear map with a nonempty output socket.
We can turn A into a multilinear form F(A) by turning
its output socket into an input socket.
Let us say that F(A) is the multilinear form of A.
We also set F(A)=A when A is a multilinear form.




 §.§ Branch decompositions and branchwidth


A branch decomposition of a (hyper)graph G consists of
(i) a tree T whose every vertex has degree either one or three, and
(ii) a bijection π between the (hyper)edge set of G and the set of
vertices of degree one in T. Deleting an edge e∈ E(T) from T
partitions T into two subtrees T_1 and T_2 that via π give rise
to a partition of the (hyper)edges of G into two sets E_1 and E_2.
The width w(e)
of the partition induced by e is the number of vertices of G that are
incident to at least one (hyper)edge in E_1 and at least one (hyper)edge
in E_2. The width of the branch decomposition (T,π) is the maximum
of the widths w(e) for e∈ E(T). The branchwidth (G) of G
is the minimum width of a branch decomposition of G.

For graphs, we recall the following upper bound on branchwidth.

[Robertson and Seymour <cit.>]
  For every n ≥ 3, (K_n) = ⌈ 2n/3 ⌉.  Consequently, (G) ≤⌈ 2 |V(G)|/3 ⌉ for every graph G.





§ EXAMPLES OF MULTILINEAR MAPS



This section reviews the specific multilinear maps that we study in this work.
Together with each map we describe its associated
tensor and a socketing of the tensor that realizes the map.



 §.§ Discrete Fourier transform




Let n≥ 2 be an integer and let ρ∈ be a primitive
n^th root of unity in the field .
[That is, ρ satisfies ρ^n=1 and ρ^n/k≠ 1
for all integer divisors k≥ 2 of n.]
Define the n× n symmetric matrix Φ with
entries Φ_i,j=ρ^(i-1)(j-1) for all i,j∈ [n].
The discrete Fourier transform (DFT) of order n at ρ is the
linear map L:^[n]→^[n] defined for all x∈^[n]
by the matrix-vector multiplication L(x)=Φ x.
In particular, the matrix Φ is the tensor
associated with L. The matrix Φ has two modes,
namely its rows and columns.
To realize L, take the columns of Φ as the input socket, and the rows
as the output socket.



 §.§ Determinant and permanent


Let us write S_n for the symmetric group of all permutations
σ:[n]→ [n]. Let c(σ) be the number of
cycles in the cycle decomposition of σ,
where each fixed point of σ is counted as a cycle.
Further standard examples of multilinear operators include the
determinant and permanent

    A    =
    ∑_σ∈ S_n(-1)^n-c(σ)∏_i∈[n]a_i,σ(i)    A    =
    ∑_σ∈ S_n∏_i∈[n]a_i,σ(i) ,

both of which are n-linear in the n columns (or the n rows) of the
matrix A=(a_ij)_i,j∈[n].
The determinant and permanent are associated with order n tensors
T̂^(), T̂^()∈^[n] × [n] ×⋯× [n]
defined for all j = (i_1,i_2,…,i_n) ∈ [n] × [n] ×⋯× [n] by


    T̂^()_j   =
    
    (-1)^n-c(j)   if  j ∈ S_n
    
    0    otherwise   T̂^()_j   =
    
    1    if  j ∈ S_n
    
    0    otherwise.

To realize the determinant and permanent using (<ref>), we socket T̂^() and T̂^() with the rows of A.
The determinant and permanent tensors are not closed under taking Kronecker products, because the Kronecker product of two tensors has the same order as its input tensors, but for each n there is exactly one determinant (permanent) tensor of order n.



 §.§ Matrix multiplication



Let n, r, and m be positive integers.
Perhaps the most fundamental example of a bilinear map is the map
that multiplies
an n× r matrix A=(A_ij)_i∈[n], j∈ [r] with
an r× m matrix B=(B_ij)_i∈[r], j∈ [m]
to obtain the n× m product matrix C=(C_ij)_i∈ [n] ,j∈ [m]
defined for all i∈[n] and j∈ [m] by

    C_ij=∑_k∈[r] A_ikB_kj .

It is natural to view the input A∈^[n]× [r] as a 2-tensor,
and similarly so for the input B∈^[r]× [m], and the
output C∈^[n]× [m]. Thus, (<ref>) is naturally
associated with the 6-tensor
T̂∈^[n]× [r]× [r]× [m]× [n]× [m] with
entries defined for all i_1,i_2∈ [n], j_1, j_2∈[m], and
k_1,k_2∈[r] by

    T̂_i_1k_1k_2j_1i_2j_2=
    
    1    if i_1=i_2 and j_1=j_2 and k_1=k_2;
    
    0    otherwise.

To realize (<ref>) using (<ref>), we can use the socketing
grouped by parentheses in
([n]× [r])×([r]× [m])×([n]× [m]), where the first
two groups are the two input sockets corresponding to A and B,
and the last group is the output socket corresponding to C.

Let us write ⟨n|,r,m|$⟩ as a shorthand for the tensor (<ref>).
From (<ref>) and the definition of the Kronecker product
it is immediate that matrix-multiplication tensors satisfy

    ⟨n|_1,r_1,m_1|⊗⟩⟨n|_2,r_2,m_2|=⟩⟨n|_1n_2,r_1r_1,m_1m_2| ⟩.

That is, matrix multiplication tensors are closed under taking  Kronecker products.




 §.§ Group algebra product


Let(𝔸,+)be an Abelian group of ordern.
Another fundamental example of a bilinear map is to convolve
two vectorsf∈^𝔸andg∈^𝔸according to the
group operation of𝔸to obtain
the vectorh=f*g∈^𝔸, defined for allk∈𝔸by

    h_k=∑_j∈𝔸f_(k-j) g_j .

The map (<ref>) is associated with the3-tensorT̂∈^𝔸×𝔸×𝔸defined for alli,j,k∈𝔸by

    T̂_ijk=
    
    1    i+j=k;
    
    0    otherwise.

A socketing of the three modes of (<ref>) that
realizes (<ref>)
is to take the first two modes as two input sockets corresponding
tofandg, respectively, and to take the last mode as the
output socket corresponding toh. The vector space^𝔸equipped with the convolution product*is the
 group algebra[𝔸].




 §.§ Kruskal operator


Letn_1,n_2,…,n_ℓandrbe positive integers.
Matrix multiplication generalizes naturally to theℓ-linear task
of multiplyingℓmatricesA^(1)∈^[n_1]×[r],A^(2)∈^[n_2]×[r],…,A^(ℓ)∈^[n_ℓ]×[r]to obtain the order-ℓtensorY∈^[n_1]×[n_2]×⋯×[n_ℓ]defined for all(i_1,i_2,…,i_ℓ)∈[n_1]×[n_2]×⋯×[n_ℓ]by

    Y_i_1i_2⋯ i_ℓ
    =∑_j∈ [r] A^(1)_i_1jA^(2)_i_2j⋯ A^(ℓ)_i_ℓ j .

This operator is known as the  Kruskal operator <cit.>
[Kolda <cit.> calls this operator the Kruskal operator. Kruskal <cit.> studied the case ℓ=3 in particular.]
of the matricesA^(1),A^(2),…,A^(ℓ).

The map (<ref>) is associated with the3ℓ-tensor

    T̂∈^
    [n_1]× [r]×
    [n_2]× [r]×⋯×
    [n_ℓ]× [r]×⋯×
    [n_1]× [n_2]×⋯× [n_ℓ]

defined for alli_1,i_1'∈[n_1],i_2,i_2'∈[n_2],…i_ℓ,i_ℓ'∈[n_ℓ]andj_1,j_2,…,j_ℓ∈[r]by

    T̂_i_1j_1i_2j_2⋯ i_ℓ j_ℓ i_1'i_2'⋯ i_ℓ'=
    
    1    if i_1=i_1', i_2=i_2', …, i_ℓ=i_ℓ' and j_1=j_2=…=j_ℓ;
    
    0    otherwise.

A socketing of (<ref>) that realizes (<ref>)
is to take each of theℓpairs of modes[n_1]×[r],[n_2]×[r],…,[n_ℓ]×[r]as an input socket and the finalℓmodes[n_1]×[n_2]×⋯×[n_ℓ]as the output socket.

Let us write⟨n|_1,n_2,…,n_ℓ|r|$⟩ for the tensor
(<ref>). Analogously to the closure property (<ref>)
for matrix multiplication tensors, we observe that

    ⟨n|_1,n_2,…,n_ℓ|r|⊗⟩⟨n|_1',n_2',…,n_ℓ'|r'|=⟩⟨n|_1n_1',n_2n_2',…,n_ℓ n_ℓ'|rr'| ⟩.

That is, Kruskal operator tensors are closed under taking Kronecker products.
Furthermore, we observe that ⟨n|,r,m|=⟩⟨n|,m|r|$⟩. That is, the
matrix multiplication tensor (<ref>) is the special case
of the Kruskal operator tensor (<ref>) whenℓ=2.




 §.§ Homomorphism-counting forms


Further multilinear operators arise
naturally by algebraization of combinatorial problems. For example,
the linear form of the matrix multiplication map

    ∑_i,j,k∈ [n]A_ijB_jkC_ki

can be used to count the number of triangles in graph, and the form

    ∑_i,j,k,ℓ∈ [n]
    A_ijB_ikC_iℓD_jkE_jℓF_kℓ

can be used to count the number of occurrences of any 4-vertex subgraph in
a graph by varying the sixn×nmatricesA,B,C,D,E,F∈^[n]×[n]. A more complicated
variant takes as input four3-tensorsA, B, C, D ∈^[n]×[n]×[n]of shapen×n×nand considers the linear form

    ∑_i,j,k,ℓ∈ [n]
    A_ijkB_ikℓC_ijℓD_jkℓ .

Ann^4-δtime algorithm for the associated trilinear map(A, B, C) ↦Dwould imply an algorithm for the Max 3-Sat problem with running time(2-δ')^n<cit.>.

The forms (<ref>), (<ref>), and (<ref>) are all special cases of what we call a homomorphism-counting form
defined by a hypergraphP, or, succinctly, aP- linear form.
In more precise terms, letPbe ak-uniform hypergraph onv ≥kvertices and write[v]kfor the set ofk-element subsets of[v].  For each hyperedgeS={i_1,i_2,…,i_k} ∈E(P)⊆[v]kofP, letX^(S)be an input tensor of shape[n]^S.  The P-linear form of  order n is the form

    ∑_σ∈ [n]^V(P)∏_S ∈ E(P) X^(S)_σ|_S,

whereσ|_Sis the restriction ofσto thekindices inS.
For example, the forms (<ref>), (<ref>), and (<ref>) are theP-linear forms corresponding toPbeing a triangle, aK_4, and a complete3-uniform hypergraph on4vertices, respectively.
It is immediate that (<ref>) is an|E(P)|-linear map.

The map (<ref>) is associated with a tensorT̂of orderk ·|E(P)|whose modes areM =
{ (S,i) | S ∈E(P), i ∈S  }with index setsJ((S,i)) = [n]. The entries ofT̂are defined by

    T̂_j=
    
    1    if ∃σ∈ [n]^V(P) such that j_(S,i)=σ_i for every (S,i) ∈ M;
    
    0    otherwise.

A socketing of (<ref>) that realizes (<ref>) is given by one input socket for each hyperedgeS ∈E(P)such that
the socket contains thekmodes(S,i)fori ∈S.

Let us observe that the tensors of the forms (<ref>) and
(<ref>) are actually the same – they are both of order12, have volumen^12, and are in fact equal to the outer product
of then ×n ×nidentity tensor with itself4times
after renaming of modes.  However, due to the difference in
socketing, the forms are computationally very different.
We show in <ref> that while there
are non-trivial tensor network algorithms for evaluating
(<ref>), no such algorithms exist for
(<ref>).

Let us write⟨n||_⟩Pfor the tensorT̂as defined
in (<ref>). Analogously to the closure
properties (<ref>) and (<ref>), we observe the
closure property

    ⟨n||_⟩P ⊗⟨n|'|_⟩P = ⟨n|n'|_⟩P .


§ TENSOR NETWORKS


 §.§ Networks


A  network (or  diagram)
consists of a finite setVof  vertices,
a finite setEof  hyperedges,
an  incidence relationI⊆V×E, and
a  boundaryB⊆E.
A network is  nondegenerate if every hyperedge is incident to at least
one vertex. In what follows we assume that all networks are nondegenerate.
A hyperedgee∈Eis a  loop ife∉Bandeis incident
to exactly one vertex.

For a vertexv∈V, let us writeI(v)={e∈E:(v,e)∈I}for the set of hyperedges incident tov. Dually, for an hyperedgee∈E, let
us writeI(e)={v∈V:(v,e)∈I}for the set of vertices incident toe.
For a networkD, we writeV(D),E(D),I(D), andB(D)to
refer to the vertices ofD, the hyperedges ofD, the incidence relation
ofD, and the boundary ofD, respectively.



    Induced networks

For a networkDand a nonempty subsetW⊆V(D),
the  induced networkD[W]consists of the vertices inWtogether with the hyperedges ofDthat are incident to at least one vertex inW; the boundary ofD[W]consists of all hyperedges that are at the boundary ofDor incident to a vertex outsideW.  Formally,

    V(D[W])   =W ,
    
      E(D[W])   ={e∈ E(D):∃ w ∈ W s.t. (w,e)∈ I(D)} ,
    
      I(D[W])   =I(D)∩(V(D[W])× E(D[W])) ,
    
      B(D[W])   =(B(D)∩ E(D[W]))∪
      {e∈ E(D[W]):∃ v ∈ V(D) ∖ W s.t. (v,e)∈ I(D)}.

For a vertexv∈V, we abbreviateD[v]=D[{v}].
Note that the boundary ofD[v]consists of all non-loop
hyperedges incident tovinD.



 §.§ Tensor networks


LetDbe a network.
We indexDby associating with each hyperedgee ∈Ean index setJ(e)of sizeℓ(e).  Induced networks inherit indexing by restriction.
Next we associate with each vertexv∈Va tensorT(v)∈^J(I(v)). We say thatDequipped with the tensors(T(v))_v ∈Vis a  tensor network.

The value of a tensor networkD, or the tensor represented byD, is a tensorT(D) ∈^J(B), defined for alli ∈J(B)by

    T(D)_i=
    ∑_j∈ J(E(D)∖ B)∏_v∈ VT(v)_ij .

Observe that in (<ref>)
the positionsiandjtogether identify a unique entry ofT(v)by projection toJ(I(v)).
We also observe that the value of a tensor network with an empty boundary is a scalar.



 §.§ Contracting tensors


LetDbe a tensor network and letW⊆V(D)be a nonempty set
of vertices. Letwbe a new element not inV.
We may  contractWinDto obtain a tensor networkD/Wby replacing the sub-networkD[W]inDwith the single vertexwwhose associated tensorT(w)is the tensor represented byD[W].  Formally,

    V(D/W)   =(V(D)∖ W)∪{w} ,
    
    E(D/W)   =E(D)∖(E(D[W])∖ B(D[W])) ,
    
    I(D/W)   =(I(D)∖ I(D[W]))∪{(w,e):e∈ B(D[W])} ,
    
    B(D/W)   =B(D) ,
    
    T(w)   =T(D[W]) .

The  cost of contractingWinDisc(D,W)=∏_e∈E(D[W])|J(e)|.
The value of a tensor network is invariant under contraction, i.e.,
for all nonemptyW⊆V(D)it holds thatT(D)=T(D/W)(see
<ref> for a proof).



 §.§ Execution and cost of a tensor network


To compute the tensorT(D)from a given tensor networkD,
we may proceed by a sequence of contractions onD. Such a process is called
executingD, and the cost ofDis the cost of a minimum-cost execution ofD. We proceed with the details.

LetD=D_0be a tensor network with at least one tensor.
Fork=1,2,…,t, select a nonempty
subsetW_k-1⊆V(D_k-1)such thatW_k-1has at least two
tensors or consists of a single tensor with a loop.
SetD_k=D_k-1/W_k-1and observe that the number of tensors and/or modes
decreases by at least one in the contraction.
Suppose thatD_tis loopless and consists of a single tensor.
We say that such a sequence of contractions is an  execution ofDintsteps.
The  cost of the execution ismax_k=1,2,…,t c(D_k-1,W_k-1).
The cost of an execution in zero steps is defined to be0.

It is immediate thatDhas at least one execution and every
execution consists of at most2|V(D)|-1steps.  By invariance under
contractions, we haveT(D_t) = T(D).
The  costc(D)ofDis the cost of a minimum-cost execution ofD.

An execution ofDof costc(D)immediately translates into an
algorithm that computesT(D)usingO(c(D)|V(D)|)arithmetic
operations in, since the contraction stepD_k =
D_k-1/W_k-1takesO(c(D_k-1, W_k-1)) ≤c(D)time to
evaluate, and there areO(V(D))steps.

LemBinaryExecutionLemma
Let D be a tensor network.
There exists a minimum-cost execution of D such that each contracted
set has size at most two. Furthermore, if D is loopless,
we can assume that each contracted set has size exactly two.


In what follows we restrict to consider looplessDonly.
Thus while a general execution may contract arbitrary vertex sets inDin each step, we may assume without loss of generality that
the minimum-cost execution has the structure of a rooted binary tree,
whose leaves are the vertices of the tensor network, and each internal
vertex is the tensor obtained by contracting its two children.



 §.§ Cost of a multilinear map


Let us now define the cost of a multilinear map via the minimum-cost
tensor networks (and socketing) for evaluating the map.
That is, the cost of a multilinear map is defined in
terms of the best tensor network that implements the map.
In more precise terms, let

    A:^J(E_1)×^J(E_2)×⋯×^J(E_ℓ)→^J(E')

be anℓ-linear map. Consider the tensorT̂(A)ofAand the
associated input socketsE_1,E_2,…,E_ℓand the output socketE'. LetD^*be an arbitrary tensor network such thatT(D^*)=T̂(A)and
the boundary satisfiesB(D^*) =E_1∪E_2∪⋯∪E_ℓ∪E'.
Modify the networkD^*as follows.
For eachk=1,2,…,ℓ, introduce a new vertex toD^*,
make the new vertex incident to each of the modes in the input socketE_k,
and associate the new vertex with a tensorX^(k)∈^J(E_k).
Remove the modesE_1∪E_2∪⋯∪E_ℓfrom the boundary ofD^*.
Let us denote the resulting network byDand call the introducedℓnew vertices the  socket vertices ofD.
We observe thatB(D)=E'andA(X^(1),X^(2),…,X^(ℓ))=T(D).
Furthermore, the costc(D)is independent of the value ofX^(k)∈^J(E_k)fork=1,2,…,ℓ.
We say thatDis a realization ofAif it can be obtained
fromAby this process, and write𝒟(A)for the set of all
tensor network realizationsDofA.

The  cost of the mapAisc(A)=min_D∈𝒟(A) c(D).
In particular, we observe that the minimum exists since the cost of a
tensor network is a nonnegative integer and the family𝒟(A)is nonempty.




§ UPPER BOUNDS ON COST


This section presents tensor-network algorithms for the maps
introduced in <ref>. We start with our key technical result
that cost is submultiplicative (Theorem <ref>,
stated formally as Theorem <ref>),
which then enables us to represent essentially the fastest known algorithms
using tensor networks, and, in the case ofP-linear forms, also
improve on earlier work as reviewed in <ref>.



 §.§ Submultiplicativity of cost


LetE_1,E_2,…,E_ℓ,E'be pairwise disjoint sets of indexed modes
such thatE_1,E_2,…,E_ℓare nonempty. Let

    A:^J(E_1)×^J(E_2)×⋯×^J(E_ℓ)→^J(E')

be anℓ-linear map.
For a positive integerk, we
define theℓ-linear mapA^⊗ksuch that its tensor satisfiesT(A^⊗k)=T(A)^⊗k.  Then

    A^⊗ k:^J(E_1)^k×^J(E_2)^k×⋯×^J(E_ℓ)^k→^J(E')^k.

Note thatT(A^⊗k) = T(A)^⊗kis thek-fold
Kronecker product ofT(A)with itself – that is, it has the same order,
but the index sets are larger – whereasT̂(A^⊗k)is
thek-fold outer product ofT̂(A)with itself – that is, its
index sets have the same sizes, but its order isktimes larger.

LetDbe a diagram that realizesAand letbe an execution
tree forD.
For each internal vertexxin(that is, a vertex obtained by contraction),
define the  amortized cost ofxby splitting into the following
three cases:

  * 
if neither of the two subtrees of x contains a socket vertex,
the amortized cost of x is 1;

  * 
if exactly one of the subtrees of x,
say, the subtree rooted at y (where x and y are adjacent in ),
contains at least one socket vertex, the amortized cost of
x is the maximum of the volume of the tensor at x and the volume of
the tensor at y;[Here, it is crucial to note that the volume of the other subtree rooted at x, only containing non-socket vertices, does not contribute directly to the amortized cost of x.]
  * 
if both of the subtrees of x contain at
least one socket vertex, the amortized cost of x is the
cost of the contraction to obtain x.


The  amortized costa()ofis the maximum of the amortized
costs of the internal vertices of. Since the amortized cost of each
internal vertex ofis at most its cost, we havea()≤c().
Furthermore, we observe that the amortized cost ofxin case (ii) above
may be strictly less than the cost of the contraction to obtainx.
In particular, in (ii) the amortized cost is defined
 not by the cost of the contraction but rather by volume.
This is because in ak^thKronecker power we can amortize
the cost of the aggregate transformation in case (ii) not with
a single contraction but with a sequence ofkcontractions.
This observation will form the heart of the proof of
Theorem <ref>.

Before proceeding with the proof, let us illustrate the
key ideas in visual terms. Let us start with the three illustrations in
(<ref>).
1.0figures/ub_tensor_rankdiag-base
Suppose the leftmost network in (<ref>)
is socketed so that the two modes at the top form the output socket,
and the four modes at the bottom form two input sockets so that modes in the
same socket are incident to the same vertex.
In the middle in (<ref>), we have adjoined two socket
vertices to the input sockets to obtain
a realizationD. On the right in (<ref>), we display an
execution treeforD. Observe that the bottom-most internal vertices
of, and the top-most internal vertex of, have type (ii). The
internal vertex in the center has type (iii). (There are no internal vertices
of type (i).) Supposing that all the modes have length at least2,
we also observe that the vertices of type (ii) have
amortized cost strictly less than their contraction cost.


Let us now consider thek^thpower
of (<ref>) visually, fork=4:

0.44figures/ub_kroneckerdiag-base-kron

The leftmost network in (<ref>) depicts thek-fold outer product of the network on the left
in (<ref>) with itself. Observe that we simply takekcopies of the
network, but that for the purposes of the visualization we have taken care to draw thekcopies of each mode together for the socketing. In the middle in (<ref>),
we have adjoined two socket vertices to the input sockets to obtain
a realizationD^⊗kofA^⊗k. On the right in (<ref>),
we display an execution treeforD^⊗k. Observe how each
of the internal vertices of type (ii) ingets expanded to a sequence
ofkinternal vertices in. This transformation fromtois the gist of the following theorem.

[Formal statement of <ref>]

Let D be an arbitrary realization of A and let  be an arbitrary
execution tree for D. For all positive integers k, we have

    c(A^⊗ k)≤ a()^k-1 c() .

Furthermore, this realization of A^⊗ k consists of
at most k · |V(D)| vertices.

Let D^* be the subnetwork of D with T(D^*)=T̂(A). That is,
D^* is the network induced by all the non-socket vertices of D.
Taking k disjoint copies of D^*, we obtain a network whose tensor is T̂(A^⊗ k). Attaching the resulting network
to tensors at sockets gives a realization of A^⊗ k.
Let us write D^⊗ k for this realization.

To establish (<ref>), it suffices to construct an execution
tree  for D^⊗ k whose cost satisfies
c()≤ a()^k-1 c(). We construct 
by rewriting  from leaves towards the root
to consider the k copies of each vertex in D^*. We start with leaf vertices which are the vertices of D^⊗ k.
We split the process into cases (i), (ii), and (iii)
as in the definition of amortized cost.
Let x be the internal vertex of  that we are currently considering.

In case (i), we perform the contraction indicated by x in each of the
k copies of D^* in D^⊗ k individually. This creates k new internal vertices
in  that are all copies of x.  We set these k vertices as the vertices
that correspond to x in the subsequent steps. Each of these contractions in  has the
same cost as the contraction indicated by x in . This cost is less or equal than c().

In case (ii), let y be the child of x in  such that the subtree rooted at y contains a socket vertex, and
let z be the other child of x in .
There is a single vertex in
 corresponding to y and k identical vertices
in 
corresponding to z. We contract these k vertices individually each with
the vertex that corresponds to y. This creates k new internal vertices
in , where we set the topmost vertex as the vertex that corresponds
to x in the subsequent steps. After the ith step, the corresponding tensor has i copies of modes of x and k-i copies of modes of y. The cost of the  contraction in the ith step is the cost of contracting y and z in  multiplied by the the volume of y to the power k-i and the volume of x to the power i-1. Since the volumes of x and y are less than or equal to a(), this cost is less than or equal to a()^k-1c().

In case (iii), let y and z be the two child vertices of x in .
By the structure of the earlier steps, we have that a single vertex in
D^⊗ k corresponds to y,
and similarly for z. We contract these two vertices. This creates one new
internal vertex in , which we set as the vertex that corresponds
to x in the subsequent steps. This tensor has k copies of modes of x. The cost of this contraction in  is the cost of the corresponding contraction in  to the k^th power, because both tensors have k copies of all modes compared to y and z. By definition, in case (iii) the amortized cost of contracting y and z is the same as the cost of contracting y and z. Hence the cost of this contraction in  is less or equal than a()^k ≤ a()^k-1 c().

This rewriting process produces an execution tree 
for D^⊗ k with c()≤ a()^k-1 c().


An immediate corollary is that tensor networks can use low rank
decompositions ofT(A)to efficiently evaluateA^⊗k.

[Submultiplicativity of low-rank executions]
  
  Let A: ^J(E_1)×^J(E_2)×⋯×^J(E_ℓ)→^J(E') be a multilinear map.  Define
  
  m =max{|J(E_1)|, |J(E_2)|, …, |J(E_ℓ)|, |J(E')|} and
  r =  T(A).  Then c(A^⊗ k) ≤max(r, m)^kmin(r, m)
  By taking a star-like network topology (as in (<ref>))
  we get an execution with a() = max(r, m) and cost c() = m ·
  r.


 §.§ Fast matrix multiplication


Let us now illustrate the use of Theorem <ref> by capturing
fast matrix multiplication with tensor networks. We start by considering square
matrix multiplication and then proceed to rectangular matrix multiplication.

We recall that the matrix multiplication exponentω(cf. <ref>) satisfies <cit.>
    ω = inf{log ⟨n|, n, n |_⟩3/log n :  n ≥ 2 } .

Above in (<ref>) we write⟨n|, n, n |_⟩3for the
3-tensor obtained from the 6-tensor⟨n|, n, n |$⟩ by flattening
the two modes corresponding the rows and columns of each of the three matrices.
That is, ⟨n|, n, n |_⟩3 has shape n^2× n^2× n^2.



For all ϵ>0 it holds that two n× n matrices may be
multiplied in O(n^ω+ϵ) operations by executing a
tensor network.


Fix an ϵ>0. By (<ref>) we can
let α,β,γ be three 3-tensors of
shape (c× c)× d for positive integer constants c≥ 2 and d
with c^2≤ d≤ c^ω+ϵ/2
such that α,β,γ decompose the matrix multiplication
tensor ⟨c|,c,c|$⟩ defined by (<ref>) as depicted below; the
indicesi_1,k_1,k_2,j_1,i_2,j_2refer to the tensor (<ref>).
The mode shared byα,β,γin (<ref>) has lengthd,
all the other modes have lengthc.
1.0figures/ub_matrixmul1ccc
For example, Strassen's decomposition <cit.> as depicted
in (<ref>) below realizes (<ref>) withc=2andd=7.
We use the numbering in magenta to indicate correspondence between
modes in (<ref>) and (<ref>).
1.0figures/ub_matrixmul2strassen-2-2-7

Let us next set up an application of Theorem <ref>.
LetA:^[c]×[c]×^[c]×[c]→^[c]×[c]be the bilinear multiplication map for twoc×cmatrices.
Observe that the tensor ofAisT̂(A)=⟨c|,c,c|$⟩.
To realize A, define two input sockets in (<ref>),
namely {i_1,k_1} and {k_2,j_1}
to obtain a realization D and an execution tree  as follows:
1.0figures/ub_tensor_rankdiag-base-2
Since c^2≤ d, the amortized cost of  satisfies
a()=d. The cost is c()=c^2d.

Let the matrices X,Y∈^[n]×[n] be given. We construct
a tensor network that multiplies X and Y.
We may assume that n≥ c. Let k be the unique integer with
c^k-1<n≤ c^k. Extend the matrices X,Y to
X,Y∈^[c^k]× [c^k] by inserting rows and columns with
zero-entries.

Since ⟨c|,c,c|^⟩⊗ k=⟨c|^k,c^k,c^k|$⟩ by
(<ref>), we have thatA^⊗kis the multiplication
map for twoc^k×c^kmatrices.
Using Theorem <ref> withDand, we obtainc(A^⊗k)≤a()^k-1c()=d^kc^2.
Moreover, the realizationD^⊗kofA^⊗kgiven by Theorem <ref> consists of|V(D^⊗k)|=O(k)vertices.
We can now associateXandYwith the
two socket vertices ofD^⊗k, taking care to
associateXwith the left socket (originating from{i_1,k_1}andα)
andYwith the right socket (originating from{k_1,j_1}andβ).
(Cf. (<ref>) for an illustration howDandin (<ref>) yieldD^⊗kand.)
ExecutingD^⊗kthen results in the product matrixA^⊗k(X,Y)=XYind^kc^2|V(D^⊗k)|≤c^k(ω+ϵ/2)+2|V(D^⊗k)|≤c^k(ω+ϵ)=O(n^ω+ϵ)operations for all large enoughn.


Let us next proceed to rectangular matrix multiplication, where our
strategy is to reduce to square matrix multiplication. Also observe that
in the caseω=2the upper bounds in the following theorem
are optimal up to the choice ofϵ>0because of the size of
the input/output.


For all ϵ>0 it holds that we may multiply an n× r matrix
with an r× n matrix by executing a tensor network in

  * O(n^ω+ϵ-1r) operations when r≥ n, and

  * O(n^2r^ω+ϵ-2) operations when r≤ n.

Fix an ϵ>0 and let α,β,γ be three 3-tensors
of shape (c× c)× d for constants c and d
as in the proof of Lemma <ref>.
Let the matrices X∈^[n]× [r] and Y∈^[r]×[n]
be given. We construct tensor networks that compute the product XY.

To establish (i), first pad X and Y using rows and columns
of zero-entries so that both n and r
become positive integer powers of c and n divides r.
This increases n and r by at most a multiplicative factor c.
We now have n=c^k and r=c^t for positive integers t≥ k.

Observe that we can compute the n× n product matrix XY by
taking the sum of r/n=c^t-k matrix products of size n× n.

Let us implement this computation with a tensor network.
Reshape
X to a (k+t)-tensor whose all modes have length c. The first k modes
index the rows, the last t modes index the columns.
Reshape
Y to a (t+k)-tensor whose all modes have length c. The first t modes
index the rows, the last k modes index the columns.

Connect X and Y into a network as displayed in (<ref>)
on the right.
1.0figures/ub_matrixmul3ccc-r-geq-n
That is,
we join t-k column modes of X with the matching t-k row modes of Y
using t-k identity matrices I_c
(to avoid degeneracy of the network if X and Y are removed).
Then we connect the remaining modes of X and Y
to the left and right sockets of a matrix multiplication network
for c^k× c^k matrices as depicted by ⟨c|^k,c^k,c^k|$⟩
in (<ref>); this matrix multiplication network is obtained
as in the proof of Lemma <ref>. The result is a multiplication
network⟨n|,r,n|$⟩ as depicted in (<ref>) on the left.
We execute this network using the structure on the right in
(<ref>) by first
contracting (with zero cost)
the identity matrices I_c with Y. We then execute the
remaining network as in the proof of Lemma <ref>.
For large enough n, the cost of the execution is at most
c^k(ω+ϵ/2)+2c^t-k+1, which translates to
O(n^ω+ϵr) operations since the network
has O(k) vertices whose contractions have nonzero cost.

To establish (ii), first pad X and Y using rows and columns
of zero-entries so that both n and r
become positive integer powers of c and r divides n.
This increases r and n by at most a multiplicative factor c.
We now have r=c^k and n=c^t for positive integers t≥ k.

Observe that we can compute the n× n product matrix XY by
taking (n/r)^2=c^2(t-k) matrix products of
size r× r.

Let us implement this computation with a tensor network.
Reshape
X to a (t+k)-tensor whose all modes have length c.
The first t modes index the rows, the last k modes index the columns.
Reshape
Y to a (k+t)-tensor whose all modes have length c. The first k modes
index the rows, the last t modes index the columns.
Connect X and Y into a network as displayed in (<ref>)
on the right.
1.0figures/ub_matrixmul4ccc-r-leq-n
That is, we join t-k row modes of X each to an identity matrix I_c
whose other mode is at the boundary. Similarly, we join matching t-k
column modes of Y each to an identity matrix I_c whose other mode
is at the boundary.
(This is to avoid degeneracy of the network if X and Y are removed.)
Then connect the remaining modes of X and Y
to the left and right sockets of a matrix multiplication network
for c^k× c^k matrices as depicted by ⟨c|^k,c^k,c^k|$⟩
in (<ref>); this matrix multiplication network is obtained
as in the proof of Lemma <ref>. The result is a multiplication
network⟨n|,r,n|$⟩ as depicted in (<ref>) on the left.
We execute this network using the structure on the right in
(<ref>) by first
contracting (with zero cost)
the identity matrices I_c with X and Y, respectively.
We then execute the remaining network as in the proof
of Lemma <ref>.
For large enough r, the cost of the execution is at most
c^k(ω+ϵ/2)+2c^2(t-k)+1, which translates to
O(n^2r^ω-1+ϵ) operations since the network
has O(k) vertices whose contractions have nonzero cost.


Let us conclude this subsection with a well-known lemma on rectangular
matrix multiplication that we can also capture with tensor networks.


  For all non-negative integers a,b,c and ϵ > 0
  it holds that we we may multiply an n^a × n^b matrix by an
  n^b× n^c matrix by a tensor network in
  O(n^max(a+b,b+c,a+c) (ω + ϵ) / 2) operations.


Let ϵ>0 be given.
By symmetry we may assume that a≤ c. Thus there are three cases
to consider, namely (i) a≤ b≤ c, (ii) a≤ c≤ b, and
(iii) b≤ a≤ c.

When a≤ b≤ c, we need to achieve O(n^(b+c)(ω+ϵ)/2)
operations. Toward this end, it suffices to multiply an n^b× n^b matrix
with an n^b× n^c matrix, which can be implemented as n^c-b
multiplications of two square matrices of size n^b× n^b.
Proceeding analogously as in <ref>, we obtain a network that
can be executed in O(n^c-bn^b(ω+ϵ/2)) operations.
We have c-b+b(ω+ϵ)≤ (b+c)(ω+ϵ)/2 as desired
since c≥ b, ω≥ 2, and ϵ≥ 0.

When a≤ c≤ b, we need to achieve O(n^(b+c)(ω+ϵ)/2)
operations. Toward this end, it suffices to multiply an n^c× n^b matrix
with an n^b× n^c matrix, and apply part (i) of <ref> to obtain
a network that can be executed in O(n^b+c(ω+ϵ-1)) operations.
We have b+c(ω+ϵ-1)≤ (b+c)(ω+ϵ)/2 as desired
since b≥ c, ω≥ 2, and ϵ≥ 0.

When b≤ a≤ c, we need to achieve O(n^(a+c)(ω+ϵ)/2)
operations. Toward this end, it suffices to multiply an n^a× n^b matrix
with an n^b× n^c matrix. An easy modification of part (ii) of
<ref> gives a network that can be executed in
O(n^a+c+b(ω+ϵ-2)) operations.
We have a+c+b(ω+ϵ-2)≤ (a+c)(ω+ϵ)/2 as desired
since b≤ a, b≤ c, ω≥ 2, and ϵ≥ 0.


 §.§ Homomorphism-counting for pattern graphs of small branchwidth


The following upper bound forP-linear forms whenPis a graph
of small branchwidth is our main result in this section.


  For any fixed pattern graph P and every ϵ > 0, there is a tensor network that evaluates the P-linear form of order n
  in O(n^2 + n^(P) (ω + ϵ) / 2) operations.

  Let (P) = w and consider any branch decomposition of P of
  width w.  Rooting this decomposition arbitrarily by subdividing
  any edge and taking the newly added vertex as root, we obtain a
  binary rooted tree T with |E(P)| leaves where the leaves are
  identified by the edges of P.  Let r be the root of T and
  for each vertex u of T, let:
  
  * C_u ⊆ V(P) (mnemonic: C for “crossing” or “cut”) be the set of
    all vertices of P that appear both in some leaf in the subtree
    rooted at u, and in some leaf outside the subtree.  By
    definition |C_u| ≤ w for all u, and C_r = ∅.

  
  * D_u ⊆ V(P) (mnemonic: D for “done”) be the set of all vertices of P that
    appear only in leaves in the subtree rooted at u, and not
    outside.  Note that C_u and D_u form a partition of the set of
    all vertices appearing in some leaf in the subtree rooted at u,
    and that D_r = V(P).

  
  * E_u ⊆ E(P) be the set of all leaves of the subtree
    of T rooted at u (recall that each leaf of T corresponds to
    an edge).  Note that E_r = E(P).

  
  * A_u ∈^[n]^C_u be the following order |C_u| tensor
    of shape n × n ×⋯× n defined for a position i ∈ [n]^C_u by
    
    (A_u)_i = ∑_j ∈ [n]^D_u∏_S ∈ E_u X^(S)_(ij)_|S,

    where ij ∈ [n]^C_u ∪ D_u is the Cartesian product of i
    and j.
  

  Note that A_r equals the value of the P-linear
  form (<ref>).
  Furthermore, for a leaf u of T corresponding to an edge {x_1, x_2} of
  P, A_u is easy to compute: if both x_1 and x_2 have degree
  at least 2, then A_u simply equals the input tensor
  X^{x_1,x_2}.  On the other hand if x_1 and/or x_2 has
  degree 1, then A_u is either a vector or a scalar, being either
  the row sums, column sums, or sum of all entries of
  X^{x_1,x_2}.  Each of these is readily computed in O(n^2)
  time by a contraction of X^{x_1,x_2} with an appropriate
  tensor of ones.

  As we shall see below, for every non-leaf vertex x of T with
  children y and z, the tensor A_x equals the contraction of
  A_y and A_z.
  Thus, the desired result would follow if we can show that the
  contraction of two siblings A_y and A_z can be computed by a
  tensor network in O(n^w (ω+ϵ)/2) operations.  We
  now proceed to establish this.




    Partition C_x into C_xy∪ C_xz∪ C_xyz where
    C_xyz = C_x ∩ C_y ∩ C_z, C_xy = C_x ∖ C_z
    and C_xz = C_x ∖ C_y (note that every element of C_x
    must appear in exactly one of C_y and C_z so these three sets
    indeed partition C_x).  Symmetrically partition C_y into
    C_xy∪ C_yz∪ C_xyz and C_z into
    C_xz∪ C_yz∪ C_xyz.  Note that C_yz⊆ D_x and that the contraction of A_y and A_z at position (i,j,k) for i ∈ [n]^C_xy, j ∈ [n]^C_xz, k ∈ [n]^C_xyz is exactly
    
    ∑_ℓ∈ [n]^C_yz A_y(i,k,ℓ) A_z(j,k,ℓ) =  A_x(i,j,k)

    Let a = |C_xy|, b = |C_yz|, and c = |C_xz|.  Split
    each mode in C_xy into two separate modes C_xy^(x) and
    C_xy^(y) where the former is used to index A_x and the
    latter is used to index A_y.  Similarly split the modes in C_xz and
    C_yz, but not the modes in C_xyz.

    The contraction of A_y and A_z can then be evaluated using
    rectangular matrix multiplication with the following tensor
    network.
    1.0branchwidth-sketchbranchwidth

    By <ref>, this network can be
    evaluated in
    O(n^|C_xyz| n^max(a+b, b+c, a+c) (ω+ϵ)/2)
    operations.  Since a + b + |C_xyz| = |C_y| ≤ w,
    b+c + |C_xyz| = |C_z| ≤ w, and
    a+c + |C_xyz| = |C_x| ≤ w, the number of operations used to
    contract A_y and A_z is bounded by
    O(n^w (ω+ϵ)/2).


 §.§ Fast Fourier transforms and fast convolution


Let us next express fast Fourier transforms as tensor networks and their
executions.
We start by observing that the Cooley–Tukey <cit.> fast
Fourier transform on_2^kcan be implemented as a tensor network.
We assumeρis a primitive(2^k)^throot of unity in
the field.


The discrete Fourier transform for the Abelian group _2^k can
be computed by executing a tensor network in O(2^kk) operations.

The case k=0 is immediate so let us assume that k≥ 1.
We construct a tensor
network whose execution multiplies a vector x∈^[2^k] with the DFT
matrix Φ in <ref> for n=2^k
to yield the result Φ x∈^[2^k].
Toward this end, let us write I_m for an
m× m identity matrix,

    H_2=[[  1  1;  1 -1 ]]

for the 2× 2 Hadamard-Walsh matrix, and
R^(k,j)∈^[2^k] for the vector obtained by concatenating
2^j copies of the vector

    [1,1,…,1_2^k-1-j,
     ρ^2^j· 0,ρ^2^j· 1,…,ρ^2^j·(2^k-1-j-1)_2^k-1-j]

for j=0,1,…,k-2.
We can (e.g. <cit.>) decompose Φ into a sequence
of 2^k× 2^k matrices

    Φ=P^(k)B^(k,k-1)T^(k,k-2)B^(k,k-2)⋯ T^(k,1)B^(k,1)T^(k,0)B^(k,0)

consisting of alternating butterfly matrices

    B^(k,j)=
    I_2⊗ I_2⊗⋯⊗ I_2_j⊗ H_2 ⊗I_2⊗ I_2⊗⋯⊗ I_2_k-j-1

and diagonal twiddle matrices

    T^(k,j)=diag(R^(k,j))

followed by a permutation matrix P^(k) that permutes the indices
in [2^k] viewed as k-bit strings by reversing the bit-order.
Since only the (j+1)^th Kronecker component in the butterfly
matrix B^(k,j) is a nonidentity matrix, and multiplication of
a vector with the diagonal twiddle matrix T^(k,j)
corresponds to pointwise (Hadamard) multiplication with the
vector R^(k,j) on the diagonal,
we observe that the sequence (<ref>)
can be represented as a tensor network D^* as depicted below (for k=5)
so that all the modes have length 2.
0.62figures/ub_fft_productdiag-phi-5
We can now connect the network (<ref>) to a data vector
x∈^[2^k] to obtain the network D below:
0.62figures/ub_fft_with_inputdiag-phi-x-5
Below we depict in red an execution tree  with cost 2^k+1
for the network (<ref>). Observe that the mode
permutation (multiplication with the matrix P^(5)) is not part of
the execution since the permutation amounts to merely rearranging
the modes.
0.62figures/ub_fft_with_execdiag-phi-x-5-exec
Since the network has O(k) tensors, the execution
(<ref>) can be carried out in O(2^kk) operations in .

The discrete Fourier transform for the elementary Abelian group _2^k can
be computed by executing a tensor network in O(2^kk) operations.

This proof is analogous to the proof of Lemma <ref> but omits the twiddle matrices
T^(k,j) and the final permutation matrix P^(k) from the
decomposition of the tensor network.


The following corollary is immediate from Lemma <ref>.


The group algebra products on [_2^k] and [_2^k]
can be computed by executing a tensor network in O(2^kk) operations whenever
2 is a unit in .

We start with [_2^k].
Let f,g∈^[2^k] be two vectors given as input. Our task is to
compute the _2^k-convolution f*g∈^[2^k].
The case k=0 is immediate so suppose that k≥ 1.
Recalling that f*g=Φ^-1((Φ f)· (Φ g)), where
“·” denotes an elementwise (Hadamard) product of two vectors of
length 2^k, let us construct a tensor network as follows.
First take the FFT of f and g using Lemma <ref>,
then multiply the resulting vectors elementwise, and
finally take the inverse FFT by
replacing ρ with ρ^-1 in Lemma <ref> and
multiplying with the diagonal matrix 1/2^kI_2^k.
Below we display (for k=5) the resulting network D that executes to f*g.
0.28figures/ub_convolutiondiag-conv
The execution of the network proceeds from right to left analogously to
(<ref>).
0.28figures/ub_convolution_with_execdiag-conv-exec
The cost of this execution is 2^k+1.
Since the network has O(k) tensors, the execution
(<ref>) can be carried out in O(2^kk) operations in .

The case of [_2^k] is analogous but replacing Lemma <ref>
with Lemma <ref> and modifying the diagrams accordingly.


 §.§ Yates' algorithm


A particularly simple use case for Theorem <ref>
occurs whenA:^[s]→^[t]is a linear map.
It is immediate that we can realizeAwith a two-vertex
network and an execution tree that has amortized costmax(s,t)and costst.

Then, Theorem <ref> immediately implies that we can evaluateA^⊗k:^[s]^k→^[t]using a tensor network
with costmax(s^k+2,t^k+2)andO(k)vertices. In particular, the
network can be executed inO(max(s^k+2,t^k+2)k)operations in.
This network in essence realizes Yates' algorithm <cit.>
for multiplying ans^k-length vector with thek^thKronecker power of ans×tmatrix to obtaint^k-length vector.

Applying the previous observation toA=H_2in (<ref>) withs=t=2,
we obtain Lemma <ref> as an immediate corollary. Similarly,
other choices of2×2matrices yield the algebraic core of currently
the fastest known algorithms for problems such as graph coloring
and its generalizations such as computing the Tutte polynomial of
a graph <cit.>.
In particular, the pair of mutually inverse2×2matrices

    Z_2=[[ 1 1; 0 1 ]] ,    
    M_2=[[  1 -1;  0  1 ]]

yield, asZ_2^⊗kandM_2^⊗k, the zeta and
Möbius tranforms for the lattice({0,1}^k,⊆,∩,∪)of all subsets of ak-element set, partially ordered by subset inclusion.
Theorem <ref> yields immediately the standard algorithms
(normally developed via Yates's algorithm)
for the fast zeta and the fast Möbius transforms via tensor networks.
These networks can then be combined as in the proof of
Lemma <ref> to yield the associated bilinear convolution maps
(multiplication maps in the semigroup algebra[({0,1}^k,⊆,∩,∪)]) to realize these maps inO(2^kk)operations. We omit the details due to similarity
with Lemma <ref>.





 §.§ Kruskal operator


We proceed to implement the Kruskal operator by Kroneckering and reduction
to fast rectangular matrix multiplication.


For all constants ϵ>0 and ℓ=1,2,…
it holds that we may evaluate the Kruskal operator
of ℓ matrices of shape n× r by executing a tensor network in

  * O(n^⌈ℓ/2⌉(ω+ϵ-1)r) operations when r≥ n^⌈ℓ/2⌉, and

  * O(n^2⌈ℓ/2⌉r^ω+ϵ-2) operations when r≤ n^⌈ℓ/2⌉.

Fix an ϵ>0 and let α,β,γ be three 3-tensors
of shape (c× c)× d for constants c and d
as in the proof of Lemma <ref>.
Let A^(1),A^(2),…,A^(ℓ)∈^[n]×[r] be given as input.
We construct a tensor network that computes the output Y of the
Kruskal operator (<ref>).

Without loss of generality we may assume that ℓ is even by
introducing a matrix A^(ℓ+1)∈^[n]×[r] filled with
1-entries and setting ℓ←ℓ+1. By inserting
rows and columns with zero-entries as necessary,
we can assume that both n and r are positive integer powers of c.
The key idea is now to Kronecker the matrices
A^(1),A^(2),…,A^(ℓ/2) and
A^(ℓ/2+1),A^(ℓ/2+2),…,A^(ℓ) in the vertical dimension only to obtain two matrices
B∈^[n^ℓ/2]× [r]
and
C∈^[n^ℓ/2]× [r],
respectively. We then multiply B and the transpose of C using fast
rectangular matrix multiplication from Lemma <ref> to obtain
Y∈^[n^ℓ/2]× [n^ℓ/2]. It is immediate from
Lemma <ref> that this
results in the operation counts claimed in (i) and (ii) as long as we
can realize the idea using tensor networks.
Kroneckering in the vertical dimension only can be realized by joining
all the  horizontal dimensions to a common mode, which becomes the
inner mode for matrix multiplication. The resulting network is
depicted below (for ℓ=6), where either (<ref>)
or (<ref>) is used for the
subnetwork indicated with ⟨n|^ℓ/2,r,n^ℓ/2|$⟩
depending on whether (i) or (ii) holds.
1.0figures/ub_kruskalkruskal-diag
In drawing (<ref>) we have made two abstractions. First,
each drawn mode in fact is a bundle of modes, each of lengthc.
Second, each mode in a bundle that is incident to one of the input matricesA^(1),A^(2),…,A^(ℓ)is in fact subdivided by inserting
an identity matrixI_cjust before the incidence to the input matrix.
The network (<ref>) is executed first by contracting
the (zero-cost) identity matricesI_c, then contractingA^(1),A^(2),…,A^(ℓ/2)andA^(ℓ/2+1),A^(ℓ/2+2),…,A^(ℓ), and finally proceeding
to execute the subnetwork⟨n|^ℓ/2,r,n^ℓ/2|$⟩ as
in Lemma <ref>.


 §.§ Clique-counting forms


For countingv-cliques, the general upper bound forP-linear forms
(<ref>) withP=K_vgives a running time ofO(n^(ω+ϵ)⌈2v/3 ⌉/2) = O(n^(ω+ ϵ)⌊v/3
  ⌋+ (ω+ϵ)2(v 3)).
In this section we give a slightly improved tensor-network algorithm,
matching the running time of Nešetřil and Poljak <cit.>.
For brevity, throughout this section we refer to theP-linear form
forP=K_vas thev2-linear form.


For all constants ϵ>0 and v=3,4,…
it holds that we may evaluate the v2-linear form of order n
by executing a tensor network in
O(n^(ω+ϵ)⌊ v/3 ⌋ + (v  3)) operations.

Fix an ϵ>0 and let α,β,γ be three 3-tensors
of shape (c× c)× d with constants c and d
as in the proof of Lemma <ref>.

Let u=v 3 and s=(v-u)/3. That is, s and u are positive
integers that satisfy v=3s+u. Let us write [a,b]=[b]∖ [a-1]
and set

    P_1   =[1,2s]∪ [3s+1,3s+u] ,
    
    P_2   =[1,s]∪ [2s+1,3s]∪ [3s+1,3s+u] ,
    
    P_3   =[s+1,3s] .

Now partition the subsets in [v]2 into three groups G_1,G_2,G_3⊆[v]2 such that
(a) G_1 consists of all subsets {x,y}∈[v]2
with {x,y}⊆ P_1,
(b) G_2 consists of all subsets {x,y}∈[v]2∖ G_1
with {x,y}⊆ P_2, and
(c) G_3 consists of all subsets {x,y}∈[v]2∖ (G_1∪ G_2). In particular, we observe that for all
{x,y}∈ G_3 we have {x,y}⊆ P_3.
Also observe that |P_1|=|P_2|=2s+u and |P_3|=2s.

Let the v2 matrices A^{x,y}∈^[n]× [n] for
{x,y}∈[v]2 be given as input. By inserting rows and
columns with zero-entries as appropriate, we may assume that n=c^k for
a positive integer k. The key idea is to start with a base construction
for an input of size c and then scale the construction up
to an input of size n using Theorem <ref>.

Let us now describe the base construction. First, we
introduce a subnetwork that forces consistency of modes for the inputs indexed
by G_1, G_2, and G_3, respectively.
In essence, for each mode
z∈ P_1, introduce a mode for each occurrence of z in G_1.
Similarly, for each z∈ P_2, introduce a mode for each occurrence
of z in G_2. Finally, for each z∈ P_3, introduce a mode for
each occurrence of z in G_3. Below we show an illustration for
v=8 and thus s=2, u=2.
0.50figures/ub_binomial_first_levv-2-subsume
In the drawing (<ref>) all modes have length c, and
each occurrence of z∈ [v] at the top is tacitly assumed to be subdivided
by an identity matrix I_c. (Otherwise the network would be degenerate.)
Next we introduce s copies of α,β,γ to force consistency
for the modes in [1,3s] between G_1,G_2,G_3. We also join individually
each of the u modes in [3s+1,3s+u] between G_1 and G_2.
Below we show an illustration for v=8. Observe from (<ref>)
and the numbering in magenta that each introduced copy
i=1,2,…,s of α,β,γ forces consistency between
the modes i,s+i,2s+i.
0.49figures/ub_binomialv-2-abc
The base network D realizing the c-uniform v2-linear form
is now as depicted in (<ref>).
We execute the network D using the following execution tree ,
where the identity matrices I_c immediately incident to the
v2 socket vertices are contracted tacitly and not drawn
in (<ref>).
0.49figures/ub_binomial_with_execv-2-exec
The amortized cost of the execution in (<ref>) is
a()=d^sc^u and the cost is c()≤ d^s+1c^u+2.
By Theorem <ref>, this translates to cost
at most a()^k-1c()=d^sk+1c^uk+2 for the
n-uniform v2-linear form. Since the network given
by Theorem <ref> has O(sk) tensors
and d≤ c^ω+ϵ/2, we conclude
that the number of operations is O(n^(ω+ϵ)s+u).
The lemma follows.


 §.§ The permanent


The following lemma observes that essentially the fastest known algorithm
for the permanent, namely Ryser's algorithm <cit.>, can be
realized as a tensor network. Here by essentially fastest known we
mean the base of the exponential running time. Sub-exponential speed-ups
to Ryser's algorithm are known, see e.g. Björklund <cit.>.


The permanent of an n× n matrix can be computed by executing
a tensor network in O(2^nn) operations.

We observe that Ryser's algorithm <cit.> for the
permanent (<ref>), namely the inclusion–exclusion expression

    A=∑_S⊆ [n](-1)^n-|S|∏_i∈ [n]∑_j∈ Sa_ij

is implementable with a star-shaped tensor network consisting of
n matrices of shape 2^n× n joined together by a common mode (of length 2^n), with the n modes of length n being the boundary of the network.
Each of the n matrices consists of the {0,1}-valued incidence vectors of
the 2^n subsets S⊆ [n], with one of the matrices arbitrarily
selected to contain signed rows determined by (-1)^n-|S|.
The input to the network consists of n vectors of length n,
namely the rows of the n× n input matrix A. The network is executed
by first executing the n matrix-vector multiplications, and then contracting
the resulting n vectors of length 2^n until the scalar A remains.


§ A LOWER BOUND FOR THE COST OF A MULTILINEAR MAP


In this section, we prove a general lower bound on the cost of
evaluating a multilinear map using tensor networks, as defined in
Section <ref>.  The lower bound is expressed in terms of the socket-width of a multilinear map, which we
now proceed to define.

LetA:^J(E_1)×^J(E_2)×⋯×^J(E_ℓ)
→^J(E')be anℓ-linear map.  A  socket-tree
ofAis a treewhoseℓ+1leaf vertices are the socketsE_1,E_2,…,E_ℓ, E'ofAand whose internal vertices all have
degree exactly3.  Associate with each edgee = {x_R, x_C}ofthe two subtrees(x_R, e)and(x_C, e)obtained by
removinge, where(x_R, e)is the subtree containingx_Rand(x_C,
e)is the subtree containingx_C.  LetL(x_R, e)be the set of leaves in(x_R, e)and letL(x_C, e)be the set of leaves in(x_C, e).

The setsL(x_R, e)andL(x_C,e)are both nonempty and together
partition the set of sockets.  Consider the flatteningM(,e)of the tensorT(A)such that the modes inL(x_R,e)index the rows and the modes inL(x_C,e)index
the columns ofM(,e).  The width ofateis the
rank ofM(, e), and the width ofisw()=max_e∈E() (M(,e)).

Let us write𝒮(A)for the set of all socket-trees of the
multilinear formA.  We define the socket-width ofAto bew(A)=min_∈𝒮(A) w().

The rest of this section is devoted to proving <ref>:

*

First, we prove that without loss of generality, we may restrict
attention to forms rather than general maps.


For any multilinear map A, it holds that c(A) ≥ c(F(A)).

We observe that A and F(A) satisfy T̂(A)=T̂(F(A)). Any network
D∈𝒟(A) can be modified to a network
D'∈𝒟(F(A)) by attaching a tensor X' ∈^J(E') to the boundary
of D. Let D∈𝒟(A) be such that
c(D)=c(A). The minimum-cost execution of D, followed
by contracting T(D) and X', is an execution of D'. Its
cost is c(A), since the cost of contracting of T(D) and X' is ∏_e ∈ B(D) |J(e)| and ∏_e ∈ B(D) |J(e)| ≤ c(A), because the last step of the minimum-cost execution of D contracted a set W with all modes e ∈ B(D) incident to W. Thus, c(A)≥ c(F(A)).


Furthermore,w(A) = w(F(A))for every multilinear mapA, sincew(A)only depends on the tensorT(A), but not on which of its coordinates (if any) is the output.
Thus it suffices to prove <ref> for
multilinear forms, which we now proceed to do.


  For any multilinear form F, it holds that c(F) ≥ w(F).

  Let D∈𝒟(F) be such that c(D)=c(F). It is a tensor
  network with empty boundary and a socket vertex S_i ∈ V(D) for
  each input socket E_i, where i=1,2,…,ℓ. Its tensor is
  T(D) = F(X^(1), X^(2),…, X^(ℓ)) where X^(i) =
  T(S_i) for i=1,2,…,ℓ.

  By Lemma <ref>, a minimum-cost execution of D
  can be represented by a rooted binary tree ,
  where the set of leaves of  are V(D) and each inner
  vertex represents the vertex obtained by contracting  its two children. Let  be the unique socket-tree of F that is obtained as a topological
  minor of . Slightly abusing the notation, we assume that the leaves of  are the socket vertices  S_1,S_2,…,S_ℓ instead of the sockets E_1,E_2,…,E_ℓ.
  To establish the lemma, it suffices to show that
   has cost at least w(), since w() ≥ w(F).

  Let e = {x_R, x_C}∈ E() be an edge of the socket tree
   with (M(, e)) = w(), and let e be an edge of
  the execution tree  in the subdivision of e appearing in
  . Without loss of generality we may assume that e is directed from the part of  corresponding to x_R towards the part corresponding to x_C (if not, simply switch names of x_R and x_C).  Define  S_R = L(x_R, e) and
  S_C = L(x_C, e).  Let W_R ⊆ V(D) be the set of
  non-socket vertices of D that appear on the same side of e in  with socket vertices S_R and let
  W_C be the set of remaining non-socket vertices of D.  See Figure <ref> for an illustration of all these definitions.
  Finally, let D' = D / S_R / S_C / W_R / W_C be the
  result of contracting each of these four sets of vertices of
  D.
  For notational convenience, we identify the four vertices of
  the new network with the four subsets S_R, S_C, W_R, W_C.

  

  Now, the tensor P = T(D'[W_R∪ S_R]) appears as an intermediate result in
  the execution ,[Note that the same is not true for the tensor T(D'[W_C ∪ S_C]).] hence the volume of P is a lower bound on the cost of .

  We group the modes of D' incident on S_R or W_R as shown in
  Figure <ref>: E_SW are all modes in D'
  incident exactly upon S_R and W_R, E_WC are all modes
  incident on W_R but not on S_R, E_SC are all modes incident
  on S_R but not W_R, and finally E_SWC are all modes incident
  upon S_R, W_R, and at least one of S_C or W_C.  Write
  E_S = E_SW∪ E_SC∪ E_SWC for the modes incident on
  S_R, and similarly E_C = E_WC∪ E_SC∪ E_SWC for
  all modes incident upon at least one of S_R/W_R and at least one
  of S_C/W_C.  Note that |J(E_C)| is precisely the volume of P
  which we aim to lower bound.

  

Define a matrix A ∈^J(E_S)×^J(E_C) as follows.  We identify its row indices i ∈ J(E_S) as being triples i = (i_SW,i_SC,i_SWC) ∈ J(E_SW) × J(E_SC) × J(E_SWC) and similarly its column indices j ∈ J(E_C) are triples j = (j_SC,j_WC,j_SWC) ∈ J(E_SC) × J(E_WC) × J(E_SWC).  Then the entries of A are
  
    A_(i_SW,i_SC,i_SWC),(j_SC,j_WC,j_SWC) =
      
        T(D'[W_R])_i_SW,j_WC,j_SWC   if  i_SC = j_SC∧ i_SWC=j_SWC,
    
        0    otherwise,

  In the case when E_S = E_SW (i.e., all modes incident on S_R connect only to W_R), A is simply a flattening of T(D'[W_R]). Recall that T(D'[S_R]) ∈∏_e ∈ E_S^J(e). Then for every j = (j_SC,j_WC,j_SWC) ∈ J(E_C), we have
  
    ∑_i ∈ J(E_S) A_i,j T(D'[S_R])_i   = ∑_i_SW∈ J(E_SW) A_(i_SW,j_SC,j_SWC),j T(D'[S_R])_i_SW,j_SC,j_SWC
       = ∑_i_SW T(D'[W_R])_i_SW,j_WC,j_SWC T(D'[S_R])_i_SW,j_SC,j_SWC
       = P_j_SC,j_WC,j_SWC = P_j

  (recall that P is the contraction of T(D'[W_R]) and T(D'[S_R])).
  Viewing T(D'[S_R]) as a row vector in ^J(E_S) we see that
  P is simply the
  vector-matrix product P = T(D'[S_R]) · A ∈^J(E_C).

  Symmetrically, for the other half of D', we can write Q =
  T(D'[W_C∪ S_C]) as a matrix-vector product Q = B · T(D'[S_C]) ∈^J(E_C) where B is a matrix corresponding to T(D'[W_S]) analogously to how A corresponds to T(D'[W_R]).

  Thus we have T(D) = T(D'[S_R]) · A ·
  B · T(D'[S_C]).  Recall that for each
  socket vertex S_i in the original network D, we have T(S_i) = X^(i).  Denoting X_R = T(D'[S_R]) and X_C = T(D'[S_C]), we get X_R = ⊗_S_i ∈ S_R X^(i) and X_C = ⊗_S_i ∈ S_C X^(i).[These identities use the fact that D is derived from a non-degenerate network D^* for T̂(F).  In particular, every mode in the network D is incident upon at least one non-socket vertex, hence all modes incident upon S_R are boundary modes in D'[S_R], and similarly for S_C.]
    Hence
  
    F(X^(1), X^(2), …, X^(ℓ)) = X_R · A · B · X_C.

  It follows that A · B is the flattening of T(F)
  to a matrix with rows indexed by the sockets in S_R and columns
  indexed by the sockets in S_C.  But this flattening is precisely
  the matrix M(, e), implying that |J(E_C)| ≥(M(, e)) =
  w(), as desired.


§ LOWER BOUNDS FOR SOCKET-WIDTH


In this section we establish lower bounds on socket-width for concrete
maps.



 §.§ Determinant and permanent


We now prove lower bounds for the socket width of the determinant and
permanent.
Let us start with the following trivial observation.


  For any socket tree  with n ≥ 2 leaves, there is an edge e = {x_R, x_C}
  such that n/3 ≤ |L(x_R, e)| < 2n/3.

  Consider all edges e = {x_R, x_C} such that |L(x_R, e)| ≥ n/3.  At
  least one such edge certainly exists since n ≥ 2. Indeed, an edge
  e = {x_R, x_C} incident to a leaf x_C has |L(x_R, e)| = n-1
  leaves.  Among these, choose an edge such that (x_R, e) is of minimal size.

  Assume for contradiction that |L(x_R, e)| ≥ 2n/3.  But then one of
  the two subtrees of x_R in (x_R, e) must have at least n/3
  leaves, and since they are smaller than (x_R, e), this contradicts
  the minimality of (x_R, e).

For every positive integer n≥ 2, the socket-width of both the determinant and the permanent of an n× n matrix is at least n⌈ n/3 ⌉.

  Let  be a socket tree for the permanent, let e = {x_R,
  x_C} be an arbitrary edge in , and let k = |L(x_R, e)| be the number of
  leaves in (x_R, e).  We now show that (M(, e)) ≥nk.

  Recall that the sockets of  are the n rows of the input
  matrix.  Without loss of generality, number the rows so that the
  leaves of (x_R, e) are rows 1 to k, and the leaves of (x_C,
  e) are rows k+1 to n.  The rows of M(, e) are indexed by
  a tuple j_R = (i_1, …, i_k) ∈ [n]^k, and the columns are indexed by a
  tuple j_C = (i_k+1, …, i_n) ∈ [n]^n-k.  The entry of
  M(, e) at position (j_R, j_C) is 1 if (i_1, i_2, …,
  i_n) is a permutation.  For any set U ∈[n]k, let u_1 <
  u_2 < … < u_k be the elements of U in ascending order and
  define j_R(U) = (u_1, …, u_k); let u_k+1 < u_k+2 < … < u_n be the elements of [n]
  ∖ U in ascending order and define j_C(U) = (u_k+1, u_k+2, …,
  u_n).

  For U_1, U_2 ∈[n]k, the entry of M(, e) at position (j_R(U_1),
  j_C(U_2)) equals 1 if U_1 = U_2, and 0 otherwise. This
  induces a nk×nk identity submatrix
  of M(, e), implying (M(, e)) ≥nk.

  By Claim <ref>,  has an edge e with n/3
  ≤ k < 2n/3.  For that edge, (M(, e)) ≥nk≥n⌈ n/3 ⌉, which completes the proof for the permanent.

  For the determinant, the only change is that some
  entries of M(, e) become -1 instead of 1, but this only changes the
  identified submatrix from an identity matrix to a diagonal matrix and
  in particular does not change its rank.


The preceding proof is similar to a lower bound by Nisan
(<cit.>, Lemma 2) used to obtain lower bounds for algebraic
branching programs.  But the lower bounds obtained there can be made
as sharp asΩ(2^n)whereas in our setting, we cannot rule out
the possibility of a tensor network that avoids splitting thenvariables in two approximately equal size parts.  This means that the best we can
obtain with our current method isn⌈n/3 ⌉instead ofnn/2.




 §.§ P-linear forms


Supposeis a socket tree for aP-linear form for ak-uniform
hypergraphPonvvertices.
Recall that the sockets of this form correspond to the elements
ofE(P) ⊆[v]k.
Given an edgee ∈E()and a vertexx
∈V(), we write 
    U(x, e)=⋃_S ∈ L(x,e) S ⊆ [v].

  Let  be a socket tree for the P-linear form of order n.
  Let e =
  {x_R, x_C}∈ E() and suppose |U(x_C, e) ∩ U(x_R, e)|= u.
  Then the socket width of  is at least w() ≥ n^u.


  Let m = n^u. We show that (M(, e)) ≥
  m by identifying an m × m identity submatrix of M(, e).

  Define E_R = { (S, i) | S ∈ L(x_R, e), i ∈ S }
  to be the modes contained in the sockets on the x_R side of e, and analogously E_C
  = { (S, i) | S ∈ L(x_C, e), i ∈ S }.

  The rows of M(, e) are indexed by J(E_R) and the columns are indexed by J(E_C). We will consider an m × m submatrix of M(, e) whose rows and columns are indexed by the elements σ∈∏_i ∈ A [n_i]. More specifically, each row of the submatrix is indexed by j_R where j_R|_(S,i)=σ|_i for i ∈ A, and j_R|_(S,i) = 1 for i ∉A. Each column is indexed by j_C where j_C|_(S,i)=σ|_i for i ∈ A and j_C|_(S,i)=1  for i ∉A.

The value of the submatrix at position (j_R, j_C) is 1 if there exists σ∈∏_i ∈ A [n_i] such that j_R|_(S,i)=σ|_i and j_C|_(S,i)=σ|_i for all i ∈ A and 0 otherwise. We obtain
  an m × m identity matrix as desired.


From this claim, the branchwidth-based lower bound is immediate.


  For any hypergraph P, the socket width of the P-linear form
  of order n is at least n^(P).

  Any socket tree for a P-linear form can be directly
  viewed as a branch decomposition of P.  Thus, by definition every
  socket tree for the form has an edge e = {x_R, x_C} where
  |U(x_C,e) ∩ U(x_R,e)| ≥(P), and the lemma now follows
  from <ref>.


For counting homomorphisms of hypergraph cliques, that is,
theP-linear form for the completek-uniform hypergraphP=[v]konvvertices, we need the
following simple lower bound on the branchwidth of complete
hypergraphs (which is most likely known, but we are not aware of a
reference).


  For v > k ≥ 3, the complete k-uniform hypergraph on v vertices has branchwidth v.

  Let T be an arbitrary branch decomposition of the hypergraph.

  First, we note that for every edge e = {x_R, x_C}, either
  U(x_R, e) = [v] or U(x_C, e) = [v].  Indeed, suppose for contradiction
  that there is an edge e = {x_R, x_C} with
  i_1 ∉U(x_R, e) and i_2 ∉U(x_C, e) for some not
  necessarily distinct i_1, i_2 ∈ [v].  Since k ≥ 2, there exists
  an edge S ⊇{i_1, i_2}, and that edge must appear in
  either (x_R, e)  or (x_C,  e).

  Now, if U(x_R, e) = U(x_C, e) = [v] for some e = {x_R, x_C},
  then we are done, so suppose for contradiction that this is not the case.
  Direct each edge towards the subtree which covers [v].  For
  example, if U(x_R, e) = [v] and U(x_C, e)  [v] then e is
  directed towards x_R.

  Since a tree is acyclic, there must be some vertex x^* such that
  all edges incident to x^* are directed towards x^*.  The vertex
  x^* cannot be a leaf, because the subtree consisting of a leaf
  contains a single edge, which covers k < v elements.  Thus x^*
  has degree 3.  Let the three edges incident to x^* be e_1 =
  {x^*, y_1}, e_2 = {x^*, y_2} and e_3 = {x^*, y_3}.  Since
  all the edges are directed towards x^* there are i_1, i_2, i_3
  ∈ [v] such that i_1 ∉S(y_1, e_1), i_2 ∉S(y_2,
  e_2) and i_3 ∉S(y_3, e_3).  Since k ≥ 3 there exists a
  edge S ⊇{i_1,i_2,i_3}.  But now the edge S cannot appear in any of the subtrees rooted at y_1, y_2, or y_3.
  Since these three subtrees together cover all leaves, this yields
  the desired contradiction and we conclude that there must exist an
  edge e such that U(x_R, e) = U(x_C, e) = [v] and therefore the branchwidth is v.


 §.§ Kruskal operator


We say that anℓ-linear Kruskal operator is n-uniform if
the lengths of the modes satisfyn=n_1=n_2=⋯=n_ℓ.


For positive integers n and r, the socket-width of an n-uniform ℓ-linear Kruskal operator is at least max(n^ℓ, n^⌈ℓ/2⌉r).

  Let  be an arbitrary socket tree for the
  operator.  One of the leaves is the output socket representing the
  tensor B of shape n × n ×⋯× n (ℓ
  times) that is obtained by applying the Kruskal operator to the input matrices
  A^(i) of shape n × r, where 1 ≤ i ≤ℓ.

  Consider the neighbor x of the output socket.  It has three
  neighbors: the output socket leaf, and two other vertices x_R and
  x_C.  Either the subtree rooted at x_R or the one rooted at x_C
  must contain u ≥⌈ℓ/2 ⌉ input sockets (since
  together they contain all ℓ input sockets).  Let e be the
  edge leading to that subtree.

  We claim that (M(, e)) ≥ n^u r.  Suppose without loss of
  generality that the input sockets A^(1), …, A^(u) are in
  the subtree rooted at x_R, and that the input sockets A^(u+1),
  …, A^(ℓ) are in the subtree rooted at x together with
  the output socket.

  Each row of M(, e) is indexed by sequences of 2u indices (i_1, j_1, …, i_u, j_u) where each
  i ∈ [n] and each j ∈ [r]. Each column is indexed by
  a sequence of ℓ + 2(n-u) indices (i_1', …, i_ℓ';i_u+1, j_u+1, …, i_ℓ, j_ℓ) where each i, i' ∈ [n]
  and each j ∈ [r].  An entry of M(, e) is 1 if
  and only if i_1' = i_1, i_2' = i_2, …, i_ℓ' = i_ℓ, and
  j_1 = j_2 = … = j_ℓ.

  For any i_1, …, i_u ∈ [n] and j ∈ [r], consider the row index (i_1, j, i_2, j, …, i_u, j) and the column
  index (i_1, i_2, …, i_u, 1, 1, …, 1; 1, j, 1, j, …,
  1, j).  The
  n^u r × n^ur submatrix of M(, e) induced by these sets of
  row and column indices is the identity matrix, thus (M(, e))
  ≥ n^u r as desired.

  For the n^ℓ lower bound, we
  instead consider the edge e' joining x with the output socket.
  Now the rows are indexed by (i_1, j_1, i_2, j_2, …, i_ℓ, j_ℓ)  and the
  columns by (i_1', …, i_ℓ').  The
  n^ℓ× n^ℓ identity submatrix is obtained by taking for
  every i_1, …, i_ℓ∈ [n], the row (i_1, 1, i_2, 1, …, i_ℓ, 1) and the column (i_1, …, i_ℓ).


§ MINIMUM-COST EXECUTION


This appendix summarizes some results on tensor network contractions and minimum-cost executions. These results are not original work, but for the sake of completeness we present them here.



 §.§ Invariance property


In the following lemma, we will show that the tensor of a network is equal to the tensor of any network that is obtained by a contraction from the original network. In particular, this implies that any execution gives the same tensor.

[Invariance]

Let D be a tensor network. For all nonempty W⊆ V(D) it holds that T(D)=T(D/W).

Let W⊆ V(D) be nonempty and let i∈ J(B(D/W))=J(B(D)).
From (<ref>), (<ref>), and (<ref>),
it follows that

    T(D/W)_i
       =
    ∑_j∈ J(E(D/W)∖ B(D/W))∏_v∈ (V(D)∖ W) ∪{w}T(v)_ij
       =
    ∑_j∈ J(E(D/W)∖ B(D/W))
    T(w)_ij∏_v∈ V(D)∖ WT(v)_ij
       =
    ∑_j∈ J(E(D/W)∖ B(D))
    T(D[W])_ij∏_v∈ V(D)∖ WT(v)_ij
       =
    ∑_j∈ J(E(D/W)∖ B(D))∑_j'∈ J(E(D[W])∖ B(D[W]))∏_w∈ WT(w)_ijj'∏_v∈ V(D)∖ WT(v)_ij
       =
    ∑_j∈ J(E(D/W)∖ B(D))∑_j'∈ J(E(D[W])∖ B(D[W]))∏_v∈ V(D)T(v)_ijj'
       =
    ∑_j”∈ J(E(D)∖ B(D))∏_v∈ V(D)T(v)_ij”
       =T(D)_i .


 §.§ The structure of a minimum-cost execution


In this section, we analyze the structure of a minimum-cost execution. In particular, we will prove Lemma <ref> that states that each contracted set has size at most two and also show that one can always contract adjacent vertices in a network.

*


If D contains loops, we may assume that
a minimum-cost execution first removes all the loops by contracting singleton
vertices incident to loops. Indeed, the cost of contracting
a singleton vertex is the volume of the tensor associated to it. Since the result of an execution is a single tensor, then every vertex has to be contained in a contracted set of an execution and none of the hyperedges incident to a vertex cannot be removed before the vertex is contained in a contracted set. Hence, the volume of any tensor in the tensor network is a lower
bound for the cost of any execution of D and we may contract singleton vertices.

So let us assume that D is loopless.
Suppose that a minimum-cost execution of D contains a contraction by
a set W={w_1,w_2,…,w_s} of size at least s≥ 3, and let w be the new vertex after this contraction. Then we can replace the contraction by W with
two contractions by W'={w_1,w_2,…,w_s-1} and W”={w,w_s} without increasing the
cost of the execution. The cost of contracting W' is less or equal than the cost of contracting W, because every hyperedge incident to W' is also incident to W. The cost of contracting W” is less or equal than the cost of contracting W, because the set of hyperedges incident to w is contained in the set of hyperedges incident to W'. We repeat this procedure until all contracted sets in the
execution have size at most two.


Two tensors in a tensor network are called  adjacent
if they are incident to a common mode.

[Execution by contracting adjacent tensors]

Let D be a loopless tensor network that is connected as a hypergraph.
Then, there exists a minimum-cost execution of D such that
each contracted set has size two and consists of adjacent vertices.

Consider a minimum-cost execution of D such that each contracted
set has size two; such an execution exists by Lemma <ref>.
If all contracted sets in the execution consist of adjacent vertices,
we are done. Otherwise, consider a contraction of two vertices, u and v,
such that u and v are not adjacent when they are
contracted to yield the vertex uv.
Consider the steps of the execution after this contraction step.
Let us call such a step a  relevant step
if it involves a descendant of the vertex uv. Let us modify the
execution as follows. First, delete the contraction of u and v from
the execution. Then, for each relevant step in execution order, replace
the descendant of the vertex uv with the current descendant of
either u or v so that the contraction becomes a contraction of
adjacent vertices whenever possible. If the descendants of u and v
become adjacent after a relevant step, contract the descendants
of u and v, and then continue the execution without further changes.
Since D is connected, the descendants of u and v
must eventually become adjacent.

We will show that this modification of the execution gives again an execution and that it has cost no larger than the original minimum-cost
execution.
In the modification of the execution, let the contraction sets containing u or a descendant of u and not containing a descendant of v be {u,w_1},{uw_1,w_2}, …, {uw_1… w_s-1,w_s}; similarly, let the contraction sets containing v or a descendant of v and not containing a descendant of u be {v,z_1}, {vz_1,z_2},…, {vz_1… z_t-1,z_t}. Here w_1,…,w_s,z_1,…,z_t can be vertices of D or vertices obtained after contraction steps.  Contracting {uw_1… w_s,vz_1 … z_t}  gives the vertex uw_1 … w_s v z_1 … z_t, which  appears also in the original execution. Hence, after certain steps the tensor networks in the original execution and in the modification are the same (after making necessary non-relevant contractions) and the modification of the original execution is also an execution.

First we consider the cost of contracting {uw_1… w_i-1,w_i} in the modified execution where 1 ≤ i ≤ s. There exists a j satisfying 1 ≤ j ≤ t such that in the original execution we contract {uw_1 … w_i-1vz_1 … z_j,w_i}. Then the cost of contracting {uw_1… w_i-1,w_i} is less or equal that the cost of contracting {uw_1 … w_i-1vz_1 … z_j,w_i}, because every hyperedge incident to the vertex uw_1… w_i-1 in the modified execution has to be incident to uw_1 … w_i-1vz_1 … z_j in the original execution. Otherwise, there would be an hyperedge in D that is incident only to vertices in {u,w_1, …, w_i-1,v,z_1, …, z_j}, and in particular both to vertices in {u,w_1, …, w_i-1} and in {v,z_1, …, z_j}. This is impossible, because then uw_1… w_i-1 and vz_1 … z_j would be adjacent vertices in the modified execution, but by assumption uw_1… w_s and vz_1 … z_t are the first adjacent descendants of u and v in the modified execution. Similarly, we can show that the cost of contracting {vz_1… z_j-1,z_j} in the modified execution for 1 ≤ j ≤ t is less than the cost of contracting a set in the original execution.

Second we consider the cost of contracting {uw_1… w_s,vz_1 … z_t} in the modified execution. Without loss of generality, we can assume that in the original execution the vertex uw_1… w_svz_1 … z_t is obtained from contracting {uw_1… w_svz_1 … z_t-1, z_t}. We will show that the cost of contracting {uw_1… w_s,vz_1 … z_t} in the modified execution is less or equal than the cost of contracting {uw_1… w_svz_1 … z_t-1, z_t} in the original execution. Indeed, every hyperedge incident to uw_1… w_s in the modified execution is incident to uw_1… w_svz_1 … z_t-1 in the original execution, because otherwise there would be an hyperedge in D that is incident only to vertices in {u,w_1, …, w_s,v,z_1, …, z_t-1}, and in particular both to vertices in {u,w_1, …, w_s} and {v,z_1, …, z_t-1}. However, by assumption uw_1… w_s and vz_1 … z_t are the first adjacent descendants of u and v in the modified execution. Similarly, every hyperedge incident to vz_1 … z_t in the modified execution is incident to uw_1… w_svz_1 … z_t-1 or z_t in the original execution, because otherwise there would be an hyperedge in D that is incident only to vertices in {u,w_1, …, w_s,v,z_1, …, z_t-1}, and in particular both to vertices in {u,w_1, …, w_s} and {v,z_1, …, z_t-1}. This contradicts that uw_1… w_s and vz_1 … z_t are the first adjacent descendants of u and v in the modified execution.

The modified execution consists of at least one less contraction of nonadjacent
vertices. Repeating this procedure until there are no contractions of nonadjacent vertices completes the lemma.


 §.§ Finding a minimum-cost execution


In this section, we will show a recurrence relation for the cost of a tensor network. To this end, letDbe a loopless connected tensor networkDwith at least two tensors.
Let us write𝒞(D)for the set of allW⊆V(D)such
thatD[W]is connected.


Let D be a loopless connected tensor network D with at least two tensors. The cost of D satisfies the recurrence

    c(D)
    =min_W∈𝒞(D)
     V(D)∖ W∈𝒞(D)max(c(D[W]),c(D[V(D)∖ W]),c(D/W/(V(D)∖ W))) .

We recall that the cost c(D) of D is the cost of minimum-cost execution of D. We can view a minimum-cost execution of D as a rooted tree
such that the root has degree two, all non-root internal vertices have
degree three, and each leaf vertex is a tensor of D.
By Lemma <ref>, the two subtrees rooted at the neighbours of the
root define two disjoint sets U,W⊆ V(D) such that U∪ W=V(D)
and both D[U] and D[W] are connected. The cost of the contraction
at the root is c(D/W/U), and the subtrees rooted at U and W
cost c(D[U]) and c(D[W]), respectively. The recurrence considers
all the partitions to connected sets, including the optimum.