Efficient Decomposition of High-Rank Tensors
Adam S. Jermyn





adamjermyn@gmail.com
Institute of Astronomy, University of Cambridge, Madingley Rise, CB3 0HA, UK


             



Tensors are a natural way to express correlations among many physical variables, but storing tensors in a computer naively requires memory which scales exponentially in the rank of the tensor.
This is not optimal, as the required memory is actually set not by the rank but by the mutual information amongst the variables in question.
Representations such as the tensor tree perform near-optimally when the tree decomposition is chosen to reflect the correlation structure in question, but making such a choice is non-trivial and good heuristics remain highly context-specific.
In this work I present two new algorithms for choosing efficient tree decompositions, independent of the physical context of the tensor.
The first is a brute-force algorithm which produces optimal decompositions up to truncation error but is generally impractical for high-rank tensors, as the number of possible choices grows exponentially in rank.
The second is a greedy algorithm, and while it is not optimal it performs extremely well in numerical experiments while having runtime which makes it practical even for tensors of very high rank.




                              
Efficient Tree Decomposition of High-Rank Tensors
    Adam S. Jermyn
    December 30, 2023
=================================================







§ INTRODUCTION


In recent years the language of tensors and tensor networks has become the standard for discussing correlations amongst physical variables <cit.>.
Disparate concepts such as lattice partition functions <cit.>, differential equations <cit.>, high-energy field theories <cit.>, and quantum computing <cit.> are now being written and examined in this language, and this trend is only increasing.

Alongside the conceptual shift towards tensors and tensor networks has come a numerical shift.
Techniques such as the density matrix renormalization group (DMRG) and matrix product states (MPS) are now recognized as special cases of tensor network manipulations and representations <cit.>.
Likewise there have been renewed efforts towards multidimensional numerical renormalization, and these lean heavily on fundamental tensor operations <cit.>.

Along the way a few key operations have been identified as bottlenecks <cit.>.
Among the more prominent are contraction of networks <cit.> and efficient representation of individual tensors <cit.>.
This work is focused on the latter.

There have been several attempts to devise efficient means of representing high-rank tensors <cit.>, and among the most successful for physical applications is the tree tensor <cit.>. Tree tensors represent a means of factoring a high-rank tensor into an acyclic network (tree) formed of rank-3 tensors <cit.>.
The advantage of the tree representation is that it may be contracted efficiently: matrix elements  of tree networks may be computed with little memory overhead.

Factoring alone does not produce a more compact representation.
To reduce the storage requirements of a tensor tree, the dimension of the bonds between constituent tensors must be restricted, either globally or using an error cutoff <cit.>.
In the latter case the network represents an approximation of the original tensor, and if the underlying correlations are weak or well-localized then the approximation may be made extremely good even with low-dimensional bonds.

The main outstanding difficulty with using tensor trees for this purpose lies in picking a decomposition <cit.>: Which variables or indices ought to go at which points on the tree?
Historically this has been answered heuristically, and said to depend on the context of the problem at hand <cit.>.
For instance, if there is a physical reason to expect two sets of variables to be weakly coupled, they can be placed on opposing sides of the tree.
More recently methods which automate the selection have been proposed <cit.> and found to perform well.
The purpose of this work is to provide two new algorithms in this vein.

In Section <ref> I review the construction of tensor trees.
I then outline a brute force algorithm in Section <ref> and a much faster greedy algorithm in Section <ref>.
Both algorithms choose decompositions based purely on the correlation properties of the tensor, with no additional physical reasoning required.
In Section <ref> I examine the results of numerical experiments and benchmarks of both algorithms.
The brute-force algorithm provides an optimal result up to truncation error.
Unfortunately it is too slow to be practical for high-rank tensors.
The greedy algorithm is practical even at very high rank while producing results which are almost as good, and so is generally much more useful.

The algorithms and experiments discussed in this work are provided in an open-source Python package.
Details of this software are given in Appendix <ref>.



§ TENSOR TREES





 §.§ Singular Value Decomposition


The fundamental building block of a tensor tree is the singular value decomposition (SVD), or equivalently the Schmidt decomposition <cit.>.
This decomposition allows us to write any matrix M of shape (n,m) in the form

    M = U Σ V^†,

where U and V are (n,k) and (k,m) unitary matrices, k=min(n,m) and Σ is a (k,k) diagonal matrix with non-negative real entries.
The diagonal entries of Σ are known as the singular values of M.
This decomposition is shown diagrammatically in Penrose notation <cit.> in Fig. <ref>.
Briefly, squares represent tensors and lines represent indices.
Where lines attached to different tensors connect those indices are to be contracted.

If M has matrix rank less than k some diagonal entries of Σ vanish, and in this case the corresponding columns in U and rows in V^† may be eliminated.
In this way the SVD allows us to produce a compressed representation of M.
Similarly, if several diagonal entries of Σ are smaller than some threshold ϵ we may eliminate them and write

    M ≈U̅Σ̅V̅^†,

where the bar indicates that these matrices have had entries eliminated.
This approximation is optimal in the sense that there is no better low rank linear approximation for a given threshold <cit.>.




As a simple example, consider the matrix

    M = 1/√(2)(
    [ 1 1; 1 1 ]).

This may be decomposed as

    U    = (1,1),
    
    V    = (1,1)
    andΣ   = (1).

If we perturb this matrix to

    M =
    1/√(2)(
    [ 1 1; 1 1 ])
    + δ/√(2)(
    [  1 -1; -1  1 ]),

then the full decomposition has

    U    = 1/√(2)(
    [  1  1;  1 -1 ]),
    
    V    = 1/√(2)(
    [  1  1;  1 -1 ])
    andΣ   = 1/√(2)(
    [ 1 0; 0 δ ]).

If δ < ϵ we may choose to truncate this decomposition by eliminating the second row and column in Σ, producing the same decomposition as that of equation (<ref>).

As a notational convenience, after truncating Σ we define

    Û≡ U √(Σ)
    andV̂≡ V √(Σ),

and likewise for the approximate matrices.
Note that the matrix square root is defined for Σ because it is diagonal.
With this notation equation (<ref>) becomes

    M = ÛV̂^†

and likewise for equation (<ref>).
Fig. <ref> then becomes Fig. <ref>.



To apply the SVD to tensors with more than two indices we transform the tensor into a matrix and decompose that instead.
More specifically, for some tensor T with indices i_1, i_2, ..., i_N we produce two groups of indices i_j_1, i_j_2, ..., i_j_M and i_j_M, i_j_M+1, ..., i_j_N-M and identify

    M_{i_j_1, i_j_2, ..., i_j_M}, {i_j_M, i_j_M+1, ..., i_j_N-M} = T_i_1, i_2, ..., i_N,

where the notation {...} means we merge (flatten) the listed indices into a single index.
This identification is depicted graphically in Fig. <ref>.
The matrix M may then be decomposed exactly according to equation (<ref>) or approximately according to equation (<ref>).
In either case after decomposition the original indices may be made separate once more, and the resulting diagram is as shown in Fig. <ref>.







 §.§ Tree Decomposition


Tree decompositions are a natural extension of the singular value decomposition, in that they simply amount to an iterated SVD.
Using the example of Fig. <ref>, a first application of the SVD produces the result in Fig. <ref>.
Each of Û and V̂^† may then be decomposed further.
For instance we may group indices i_1 and i_2 together so that

    M_{i_1, i_2}, {i_3, i_i} = Û_i_1, i_2, i_3, i_i,

where i_i is the internal line.
The matrix M may then be decomposed with the SVD.
The result is shown in Fig. <ref> where A and B result from the decomposition of U.

This process may be iterated until there are no tensors in the graph with more than three indices.
This choice was made because the cost of storing tensors scales exponentially in their rank and three is the minimal rank which allows a tree structure[In some cases this may incur unnecessary overhead. For instance if the indices of a tensor of rank four are maximally correlated attempting to decompose the tensor into two tensors of rank three doubles the storage required. However this cost is a constant overhead and so is typically acceptable.]
Furthermore the choice to proceed until all tensors have three indices is the only choice which can represent tensors of any rank greater than two without invoking tensors of varying ranks, and this uniform structure is considerably simpler to manipulate than a more general one with variable tensor ranks.

At this point further decomposition may be warranted on the original indices, but further decomposition on indices which are internally contracted just introduces further intermediate tensors with no corresponding reduction in rank.
This is because the SVD is already optimal, and so if an index results from performing the SVD it cannot be further compressed without raising the error threshold.
Because the SVD just introduces intermediate tensors between sets of indices it cannot introduce cycles into the graph and so the result of this iterative decomposition is a tree.
What we mean by a tensor tree then is an acyclic tensor network (graph) composed of only tensors with three or fewer indices.





 §.§ Choice Sensitivity


To understand why the tree decomposition is sensitive to the precise choice of tree, consider the tensor

    T_ijkl = e^ije^kl.

There are three possible tree decompositions of T, shown in Fig. <ref>.
The corresponding singular values for each are plotted in Fig. <ref>.
Because T factors precisely into a product of terms dependent on i and j and ones dependent on k and l, the first decomposition has only one non-zero singular value.
In the numerical decomposition the second largest singular value is a factor of 10^16 smaller than the largest, which is within the floating point tolerance.
By contrast the other two decompositions have many non-zero singular values, and may require several of these to be included to produce a good approximation.
So for instance an approximation which is good to one part in 10^4 requires five singular values to be retained with these decompositions.

Picking a good decomposition only becomes more important for tensors with more indices.
This is because each internal index must encode the correlations between all indices on either side of it.
Separating a highly correlated pair of indices therefore requires increasing the rank of every intermediate index between them to accommodate their correlations.







§ ALGORITHMS




 §.§ Brute Force



A straightforward way to optimize over a finite set of choices is to evaluate every choice and pick the best one.
For tensors with only a few indices this is possible: we simply enumerate all possible tree decompositions, evaluate each one, and pick the one which requires the storing the fewest elements for a given error threshold ϵ.
This enumeration may be done quite efficiently <cit.>, but unfortunately the number of trees grows rapidly with the number of indices.
The number of trees with N leaves and three links on each internal node is

    S_N = (2n-5)!! = 2^2-N(2N-4)!/(N-2)!

<cit.>.
For large N,

    S_N∼ 2^-N(2N)!/N!.

Applying Stirling's approximation we find

    S_2N-2,N∼ 2^-N(2N/e)^2N(N/e)^-N = (2N/e)^N

which grows exponentially with N.
As such this method is only practical for tensors with a small number of indices.

Our implementation of this algorithm iterates over all trees.
It identifies and performs a sequence of SVD operations to produce each tree.
When the truncation error is negligible the order of SVD operations does not matter.
In testing however there were rare cases in which the order of application of the SVD mattered.
This may occur, for instance, when a singular value lies very near to the threshold and so can be pushed to one side or the other by truncation error.
Still, in almost all cases the brute force algorithm identifies the optimal tree decomposition, and so we do not further examine all possible SVD orderings.
As it stands the algorithm is practical only for relatively small tensors, and examining SVD orderings would only make this worse.



 §.§ Greedy Choice



While considering every possible tree decomposition ensures that we find the best one, there are useful heuristics which can often do nearly as well.
One such heuristic is to construct the tree by locally minimizing the number of terms left in the truncated SVD at each stage.
More specifically, we consider all factorizations of the tensor T with indices i_1,i_2,...,i_N of the form

    T_i_1,i_2,...,i_N≈ A_i_j_1,i_j_2,i_i T'_i_i,i_j_3,...,i_j_N,

where there is summation over the repeated internal index i_i.
Thus we are looking for the pair of indices to factor away from the rest which minimize the rank of the resulting internal index.
This requires checking N(N-1)/2 factorizations except in the special case where N = 4, where symmetry allows us to check just 3 cases.
Once we have found the best pair to factor away from T we may iterate and search for the best index pair to factor out from T'.
This process ends when the process returns two tensors each with three indices.

It is possible to further improve the performance of this algorithm by noting that most indices of T' are the same as those of T.
If T is exactly represented by the contraction of A and T' then we do not need to check pairs of these again, as they represent precisely the same decomposition as before.
If T is only approximately represented then in principle the rank of these pairs may have changed, but in practice for a sufficiently small error threshold ϵ we expect this to be unlikely.
As a result we only need to examine pairs involving the new index i_i, of which there are N-2.
Iterating until there are just three indices on each tensor requires

    1/2N(N-1) + ∑_j=1^N-4 (j+2) = N^2 - 2N - 2

pairs to be tested if N > 4.
In the limit of large N this scales quadratically, and so is immensely preferable to the exponential scaling of the brute force algorithm.

Note that for simplicity our implementation of this algorithm does not include this improvement, and has cubic runtime in N.
Regardless the number of comparisons is polynomial in N rather than exponential, which is the key advantage of this algorithm.

It is important to emphasize that this algorithm does not generally produce the optimal tree decomposition, though in specific instances it may.
For instance if the tensor was formed from an outer product of rank-2 objects this algorithm would correctly identify the pairs of indices which are uncorrelated from all others and the resulting tensor tree would be composed of the original rank-2 objects, augmented by links of dimension one to one another and possibly rescaled.
Likewise if the tensor was formed by contracting a sequence of identical rank-3 tensors against one another, as in the N-point correlation function of the one-dimensional Ising model (see Section <ref>), the entanglement between any pair of indices and the remainder of the indices is guaranteed to be monotonic in the distance within the pair and hence a greedy approach correctly identifies that adjacent indices ought to be split off from the rest.
These scenarios are quite special, however, and do not necessarily generalize.
Hence the greedy algorithm should be treated as producing approximations to the best choice of decomposition rather than the best decomposition outright.



 §.§ Runtime Comparison


In both algorithms the runtime of the SVD typically scales exponentially in N.
This is because the number of entries in the matrix scales as d^N, where d is the typical dimension of the indices.
This may be avoided by employing approximate decompositions and rank estimations <cit.>.
Nevertheless such approximations are not our focus and so we proceed with the naive exponential implementation.

Because the runtime is proportional to the cost of performing decompositions and/or rank estimations, for comparing these two algorithms we may just examine how many such operations must be performed.
In the brute-force algorithm this scales as (2N/e)^N, while in the greedy algorithm it scales as N^2.
For N > 1 the latter is smaller, and even for modest N the difference in runtime may be enormous.
For example when N=10 the brute force algorithm requires approximately (20/e)^10≈ 5× 10^8 comparisons while the greedy one just requires of order 10^2 or 10^3 comparisons depending on the implementation.
This makes the latter preferable so long as it produces results of comparable quality.



§ NUMERICAL EXPERIMENTS



We have already discussed the relative runtime of our algorithms and so we proceed directly to examining their relative performance on a variety of tasks.
For simplicity our implementation of the greedy algorithm does not take advantage of the fact that index pairs repeat between iterations.
Furthermore we take the output from the brute force algorithm to be the standard of comparison for the greedy algorithm, as it is optimal in all but certain rare cases [Because the brute force algorithm relies on the iterated application of SVD it may fail to produce optimal results if approximation error from an early SVD application suffices to change the effective ranks of subsequent applications. In such an instance this error could result in a sub-optimal tree being selected. With ϵ set as small as in our numerical experiments this ought to be rare, and indeed we only see it occurring twice. Both instances appear as the greedy algorithm outperforming the brute force one in Figure <ref>.].
Timing was done on an Intel IvyBridge processor using the Intel MKL linear algebra backend.
The linear algebra was parallelized across cores when this provided a speed improvement (i.e. for large tensors) but the search algorithms were run serially.
This prevents competition between threads for linear algebra resources.

For the threshold we use

    ϵ^2 = 10^-6∑_i Σ_ii^2.

This ensures that the relative L_2 error incurred with each eliminated singular value is no more than 10^-6, and means that the normalization of the tensor does not impact our choice of decomposition.

To begin we consider the example of equation (<ref>).
We tested the case where all indices range from 0 to 3 inclusive and found that the greedy algorithm identified the optimal decomposition.
This is not surprising, as the greedy algorithm reduces to the brute force algorithm when there are just four indices.

We further considered the case

    T_ijklm = e^ije^kle^m,

where once more all indices range from 0 to 3 inclusive.
Once more the greedy algorithm reproduces the result from the brute force algorithm, which is shown in Figure <ref>.
Indeed the links between A and B and between B and C are one-dimensional in this example, reflecting the total factorization of equation (<ref>).
This results in an overall compression from 4^5=1024 entries to just 36, and is an example of the best case for this sort of tensor compression.

The decomposition in Figure <ref> remains optimal even for more complex tensors like

    T_ijklm = e^ije^kle^m + sin(ij)cos(kl)tanh(m)

because the required bond dimension to encode the correlations between m and the remaining indices is less than that required to encode those between either of i and j or k and l.
The greedy algorithm once more performs optimally in this case.



To construct more challenging examples we next considered the family of tensors with N indices of dimension 2 where each entry is randomly and independently chosen as either 1 with probability p or 0 with probability 1-p.
In the low-p limit tensors drawn from this distribution exhibit non-uniform correlations across all indices because there are few correlations overall owing to most entries vanishing, and so clusters and hence correlations are typical.
Figure <ref> shows the factor by which each algorithm was able to compress tensors generated in this fashion as a function of the number of indices N as well as the time required to find this decomposition.
For the cases where both algorithms have tractable runtimes they perform very similarly, with the greedy algorithm usually finding optimal or near-optimal solutions.
For larger N the brute force algorithm becomes intractable but the greedy algorithm remains performant.
This can be seen in the timing curves, where the difference begins at several orders of magnitude and rises rapidly in N even on a logarithmic scale.
Beyond N=8 it was not practical to evaluate this algorithm with enough samples for robust averaging, and so the brute force benchmark ends there.



One feature of Figure <ref> worth noting is that this class of tensors becomes less compressible as N increases, and for some N the decompositions are actually less efficient than just directly storing the original tensor.
This is not just a feature of the greedy algorithm, as the optimal result shows the same trend as far as it has been evaluated.
To understand this note that for a given pair of indices the tensor may be viewed as a collection of matrices, each of which possesses those two indices.
In the large-N limit the number of such matrices increases rapidly.
This increase in the mutual entropy between any pair of indices increases the SVD rank required to maintain the specified accuracy.
When this rank approaches the original matrix rank the SVD can actually increase the number of entries being stored from nm up to (n+m)min(n,m), where n and m are the matrix dimensions.
This results in an overall increase in the memory required to store the tensor.

This incompressibility is actually a common feature of random tensors.
For instance Figure <ref> shows the same plot as in the top panel of Figure <ref> but for tensors whose entries are drawn randomly and independently from a normal distribution.
For all N the compression ratio is quite poor, in stark contrast to the more structured examples.
The two algorithms perform identically for N < 7, and for N=8,9 the brute force algorithm is somewhat better.
At larger N the compression ratio of the greedy algorithm improves and plateaus in performance near a compression ratio of unity.
We expect that unbiased random tensors ought not to be compressible so this result is likely close to optimal.



A case with more physical content is that of a tensor encoding the N-point correlation function for a lattice system.
For instance in the 1D classical Ising model with open ends states are weighted as

    P({s_i}) ∝ e^-J ∑_i=0^N-1 s_i s_i+1,

where each s_i is drawn from {0,1}.
Interpreting each s_i as an index produces a tensor with N indices, shown in Figure <ref>.
This object ought to be quite compressible because the correlations it encodes are local and so factor more readily.





Figure <ref> shows the tree decomposition produced by the brute force algorithm for the Ising N-point function with N=9 and J=1.
Note that the structure is regular and local, with indices which are neighbours in the Ising model appearing near one another in the tensor tree.
This reflects the local nature of correlations in the Ising model.

Figure <ref> shows the compression ratio of the Ising N-point function as a function of system size.
The two algorithms produce identical results, so in these cases the greedy algorithm finds the optimal tree decomposition.
Once more the brute force algorithm proved impossible to benchmark beyond moderate N and so only the greedy algorithm is shown for larger values.

Interestingly, the compression ratio in Figure <ref> decreases rapidly as a function of N, in keeping with the intuition that Ising correlations are highly redundant.
In particular, the decrease is exponential in N.
To understand this note that the tensor size is 2^N.
For finite J there is a corresponding correlation length ξ(J).
The size of the compressed tensor is set by the dimensions of the internal bonds in the tree.
This in turn is set by the number of degrees of freedom above the given error threshold, which is proportional to the entanglement entropy between the two sides of the bond, and so scales as 2^ξ(J)[For a more detailed version of this argument see e.g. <cit.>.]
As such the compression ratio scales like 2^ξ(J) - N, which is consistent with the data in Figure <ref>.



Figure <ref> shows the same compression ratio and timing data as a function of the site-site coupling J.
The timing is largely independent of J, which indicates that it is mainly being set by the combinatorics of the problem and not by the underlying matrix transformations.
This highlights the need for efficient means of choosing a decomposition.

Once more the two algorithms produce identical results for the compression ratio, with variation only near the transition in compression ratio at J=-2.
The increase in compression ratio with J occurs because as J increases correlations become longer-ranged, and so become more difficult to encode in a local model such as a tree decomposition.
This manifests as a sharp transition because near J=-2 a specific singular value crosses the threshold ϵ and becomes relevant.






As a further test consider the 2D Ising probability function

    P({s_i}) ∝ e^-J ∑_⟨ ij ⟩ s_i s_j,

where ⟨ ij ⟩ indicates all nearest-neighbour pairs.
Figure <ref> shows the compression ratio and timing for this function, again interpreted as a tensor indexed by {s_i}, as a function of J for the case of a 3× 3 open-boundary square lattice.
In addition to the brute force and greedy algorithms a variant of the greedy algorithm which uses variable-size tensors was included.
This variant leaves tensors of rank greater than three in the final decomposition if the next step of their decomposition would not produce net compression.
In that sense it is also greedy in the space of possible final tensor ranks.

Once more the algorithms produce very similar results.
In this case though the compression ratio is worst for intermediate J.
This is matched with an increase in computation time, which makes sense because the SVD in this region is more time consuming.
At extreme J the complexity falls because as the J→±∞ the ferromagnetic and antiferromagnetic ground states come to dominate the probability distribution.
This makes the tensors more compressible.

For most J the regular greedy algorithm outperformed the variable-size one.
This reflects the fact that the latter may determine that the next decomposition step produces no net improvement, or indeed produces a worse compression ratio, than doing nothing at all.
For random tensors, or indeed for the tensor in equation (<ref>) near J=-1 and J=5, this is a correct determination, but in other cases it can leads to premature termination of the decomposition and worse overall performance.

It is worth noting that in this case both greedy algorithms actually outperform the brute force algorithm in small regions near J=-1, J=0 and J=7.
These are among the cases where the order of the SVD matters because of the truncation error.

A key takeaway from these examples is that the structure of the tensor matters.
For instance it is not surprising that models with local structure such as these should be good use cases for the greedy algorithm because it is designed to take advantage of the existence a hierarchy of correlations, identifying the most tightly-correlated indices and moving outward from there.

Finally, for comparison purposes consider Example 6.4 of <cit.>.
Let

    g(x,y) ≡1/|y-α x| + 1

for fixed α.
Furthermore given N∈ℕ let

    ξ_i_1,...,i_N≡∑_μ=1^N 2^μ-N i_μ - 1,

where each i_k ∈{0,1}.
With this they define the tensor

    A_i_1,...,i_2N≡ g(ξ_i_1,...,i_N, ξ_i_N+1,...,i_2N).

This example is non-trivial because neither the structure of g nor the bitwise indexing of ξ permit any obvious factoring.

Table <ref> shows the effective dimension defined in <cit.> for both their adaptive agglomeration scheme and the greedy one introduced in this work.
To match the parameters used in their tests the dimension was set to N=8, the tolerance ϵ = 10^-8, and the test was performed for each α∈{0,0.25,0.5,0.75,1.0}.
Furthermore to account for possible differences in how the methods determine the SVD tolerance the greedy algorithm was also run with a tolerance of 10^-9.

In three cases, namely α = 0, α=0.25 and α=1.0 the greedy algorithm outperformed the agglomeration algorithm.
This was true both for ϵ=10^-8 and ϵ=10^-9.
In the remaining cases the agglomeration algorithm outperformed the greedy one.
In no case was the performance of one substantially different from that of the other.
Note also that the runtime complexities of these methods both require a number of SVD or other decompositions which scales as N^2.
Hence in both runtime and performance, at least on this test, these methods are comparable.





§ CONCLUSIONS


High-rank tensors are increasingly common in physical computation, and so there is now a demand for algorithms which can manipulate them efficiently.
I have introduced and tested two new algorithms for identifying good tree decompositions of tensors.
The first is a straightforward brute-force approach which, up to truncation error, to produces the optimal result at the expense of an exponential number of trials.
The second is a greedy algorithm which, while not guaranteed to produce the optimal outcome, performs very well in all cases examined while only requiring a quadratic number of comparisons.
These fill a crucial niche among the building blocks for tensor computations.



§ ACKNOWLEDGEMENTS


I am grateful to the UK Marshall Commission for financial support as well as to Milo Lin and Ravishankar Sundararaman for helpful discussions regarding this work.
This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.





§ SOFTWARE DETAILS



The software used in this work is available under a GPLv3 license at <github.com/adamjermyn/TensorDecomp>, and consists of a collection of Python scripts.
For simplicity the implementation of the greedy algorithm does not currently take advantage of the possibility of reusing comparisons.
The brute force algorithm enumerates all possible trees by first constructing a tree with one internal node and three indices and then iteratively inserting nodes attached to additional indices in place of edges.
These were run on Python v3.6 with NumPy v1.12.1 <cit.> and the plots were created with Matplotlib v1.5.3 <cit.> and the ggplot style, though no specific features of these releases were used.



