Wolfram Language Paclet Repository
Community-contributed installable additions to the Wolfram Language
Contraction Paths and Execution
Once a tensor network is built, two questions remain before a number falls out: in what order to merge tensors, and what numerical engine to drive each merge. The paclet separates those concerns. A is a pure ordering, found by or ; an execution engine – selected by the Method option to – turns each ordered pair into a numerical operation. This tutorial covers the path representations, both planners and their objectives, the five execution back-ends, and the symbolic Inactive form that lets you inspect a plan before evaluating it.
In[155]:=
Needs["Wolfram`TensorNetworks`"]
1. What a Contraction Path Is
A contraction path is a list of integer pairs {{i,j}, …}. Each step removes the tensors at positions i and j from the working list and appends their pairwise contraction – the opt_einsum convention, where both source positions disappear and the result lands at the end. After n-1 steps an n-tensor network collapses to a single tensor. The path is independent of the tensors themselves: the same path can drive the contraction of any network with the same hyperedge structure, which is precisely what makes a separate planner-then-engine pipeline worthwhile. To make this concrete, start with a three-tensor matrix-product-state-like network: a rank-2 left boundary, a rank-2 middle (so it stays in two dimensions), and a rank-2 right boundary, sharing two bond indices. Three tensors is the smallest case where path choice matters – for two tensors there's only one merge, and the planner is trivially correct.
In[156]:=
mps=
[{ArraySymbol["A",{2,3}],ArraySymbol["B",{3,3}],ArraySymbol["C",{3,2}]},{{1,2},{2,3},{3,4}}]
 | TensorNetwork |
Out[156]=
TensorNetwork
 |  |
|
Indices 2 and 3 are the contracted bonds; indices 1 and 4 are the free legs. The hyperedge lists {1,2}, {2,3}, {3,4} tell the paclet, per tensor, which integer label sits on each leg – no separate index-naming pass is required. Two natural pair-paths exist for this chain: {{1,2},{1,2}} contracts A with B first, then the AB result with C, and {{2,3},{1,2}} contracts B with C first, then with A. For a length-3 chain both cost the same; for longer chains the optimal answer is no longer obvious, and that's where the planner functions earn their keep. checks that a list has the right shape – every step is a list of one or two integers – and checks the stricter form in which every step is a pair (no unary identity steps from a binarization pass):
In[157]:=
path={{1,2},{1,2}}
Out[157]=
{{1,2},{1,2}}
In[158]:=
[path],
[path]
 | PathQ |
| CanonicalPathQ |
Out[158]=
{True,True}
The first {1,2} merges positions 1 and 2 of the working list {A, B, C} to yield {C, AB} – position 3 slides into position 1 (only one element was left in front), and the merged AB tensor is appended at the new end. The second {1,2} then merges C with AB to a single tensor. The bookkeeping is mechanical; both planners produce paths in exactly this form. One consequence worth internalizing early: {1,2} in step k refers to positions in the working list after the k-1 prior steps have re-indexed it – not to positions in the original tensor list. The opt_einsum convention keeps the bookkeeping local but means paths cannot be reordered or sliced without renumbering. The convenience function performs exactly this renumbering, which is why the tree form (a section below) is sometimes easier to manipulate.
2. The Greedy Heuristic
walks the network one pair at a time, always picking the locally cheapest merge. "Cheap" defaults to the size of the resulting tensor, optionally penalized by a memory term. At every step the algorithm considers each pair of tensors that share at least one index, scores the merge by the size of the resulting intermediate, and commits to the lowest-scoring pair before moving on. The cost is linear in the number of candidate pairs per step, so a network of N tensors runs in roughly O(N^2) total – fast enough to use as the default once a network grows beyond about 20 tensors, where the exponential dynamic-programming search of becomes infeasible. The greedy planner forgoes any guarantee of global optimality, but on most tensor networks the locally cheapest choice happens to be close to optimal, and the gap shrinks further with the simulated-annealing options below. Build a deterministic five-tensor test network from a random graph to see the planner at work:
In[159]:=
tn=BlockRandomSeedRandom[42];
[RandomGraph[{5,8}],3]
 | RandomTensorNetwork |
Out[159]=
TensorNetwork
|  |
|
The result is a connected five-tensor random network with mixed ranks – some legs are shared bonds and two are free legs that will survive the contraction. The first call asks for a path with all defaults; the algorithm is deterministic when no temperature is supplied, so repeated calls return the same path:
In[160]:=
 | GreedyContractionPath |
Out[160]=
{{1,2},{3,4},{2,3},{1,2}}
Four options shape the search. "Temperature" perturbs the greedy choice with simulated-annealing noise (zero or is deterministic; positive values widen the search by sometimes accepting a suboptimal pair, which lets the planner escape local minima at the cost of reproducibility). "RandomSeed" makes a noisy run reproducible by fixing the underlying random stream. "MaxNeighbors" caps how many candidate pairs are scored at each step, trading quality for speed on very dense graphs – useful for thousand-tensor PEPS-like networks where the per-step candidate list is itself large. "MemoryWeight" adds a penalty proportional to peak intermediate memory, which can change the chosen path when the cheapest-size step would build a large transient tensor. A practical pattern is to run the planner three times – with "Temperature" 0 (default), with a small positive temperature and a fixed seed, and with "MemoryWeight" 2.0 – then pick whichever path has the lowest objective on the actual network:
In[161]:=
| GreedyContractionPath |
Out[161]=
{{4,5},{2,4},{2,3},{1,2}}
The memory penalty has rerouted the planner: instead of starting with positions 1 and 2, it now begins at 4 and 5 and keeps the working tensors smaller through the run. On the small test network the difference in cost is modest; on a thousand-tensor lattice it can be the difference between contractible and out-of-memory. A reasonable rule of thumb: when peak memory matters more than wall-clock time – often the case on a laptop or for very tall MERA stacks – set "MemoryWeight" 1.0 or above. When tensors are uniformly small but there are many of them, leave the default and accept that two adjacent tensors might briefly fuse into something fat. Greedy paths are also stable input for : none of the downstream execution machinery cares whether a path came from the greedy heuristic or the optimal planner, only that it satisfies .
The "MaxNeighbors" option caps the size of the greedy beam – the number of candidate contractions considered at each step. Lower values are faster but less optimal.
In[162]:=
| GreedyContractionPath |
Out[162]=
{{1,2},{3,4},{2,3},{1,2}}
The "PreSimplify" option (default ) preprocesses the network – fuses duplicate indices and removes trivial dimensions – before the greedy search.
In[163]:=
| GreedyContractionPath |
Out[163]=
{{1,2},{3,4},{2,3},{1,2}}
Different "RandomSeed" values explore different greedy paths; with non-zero "Temperature" this is a stochastic local search.
In[164]:=
| GreedyContractionPath |
Out[164]=
{{1,2},{3,4},{2,3},{1,2}}
In[165]:=
| GreedyContractionPath |
Out[165]=
{{1,2},{3,4},{2,3},{1,2}}
Greedy also accepts a , an association of "Dimensions" / "Indices" / "Contractions" data, or an inactive TensorContract[TensorProduct[…], …] expression – useful when the caller already has a partial pipeline assembled.
3. Dynamic-Programming Optimum
runs a dynamic-programming search over subsets of tensors and returns a provably optimal path for whichever objective the Method option names. Conceptually, the algorithm enumerates subsets of the working list and asks, for each subset, what its best contraction order is in terms of the subsets it can be split into; the answer for the full set then composes the optimal answer for the whole network. The search is exponential in the number of tensors – roughly O(3^N) without pruning – so it is appropriate for networks up to about twenty tensors, beyond which the greedy heuristic of the previous section is the practical choice. The default objective is "size", which minimizes the largest intermediate tensor across the entire contraction:
In[166]:=
| OptimalContractionPath |
Out[166]=
{{1,2},{1,2},{2,3},{1,2}}
Six Method values are available, each minimizing a different scalar functional of the path. "size" (above) bounds the peak intermediate – the implicit proxy for memory in the original opt_einsum library. "flops" minimizes the total floating-point operation count summed over all pairwise merges. "max" bounds the size of the largest intermediate as a tighter peak-memory proxy than "size" (which can also be sensitive to how rank-1 boundaries are handled). "write" minimizes the total bytes written across all intermediates, which approximates cache traffic better than raw FLOPs on bandwidth-bound hardware. "combo" is a weighted blend of size and flops, useful when neither objective dominates in isolation. "limit" minimizes the per-step largest tensor under a hard budget constraint and is the right choice when an out-of-memory error must be avoided at any cost. The size and flops paths differ on this network:
In[167]:=
| OptimalContractionPath |
Out[167]=
{{1,2},{3,4},{2,3},{1,2}}
Step two flips from {1,2} (size-optimal) to {3,4} (flops-optimal): the size objective accepts a wider intermediate that costs more arithmetic, while the flops objective spreads the work across two independent left-and-right partial products before joining them. Which is preferable depends on the bottleneck – memory or compute. For modern dense-tensor contractions on a workstation with abundant RAM, "flops" usually wins on wall-clock time; for memory-constrained machines or for sparse intermediates, "size" or "limit" is safer. For very large search spaces "PruningThreshold" caps the dynamic-programming table by discarding partial solutions whose cost already exceeds the threshold, trading optimality guarantees for tractability. With pruning, the planner stops being provably optimal but still returns a path – one that is at least as good as the threshold, and often significantly better:
In[168]:=
| OptimalContractionPath |
Out[168]=
{{1,2},{1,2},{2,3},{1,2}}
On this five-tensor test the pruned search recovers the same path as the unpruned size-optimal search, because the threshold is large enough not to discard any winning branch. On a twenty-tensor network the pruning option is what makes the dynamic-programming planner usable at all; without it the table can hold tens of millions of entries. The cost-cap value is in the same units as whichever Method objective is in use, so a threshold of 100.0 against "size" means "discard any partial solution whose largest intermediate already exceeds 100 elements". Set the threshold generously at first, then tighten it once you see how the planner responds.
The Method option supports six cost functions:
Method | What it minimizes |
"size" (default) | the size of the largest intermediate tensor |
"flops" | the total floating-point operation count |
"max" | the size of the single largest tensor created |
"write" | the total write cost (sum of all intermediates) |
"combo" | a combined size+flops cost |
"limit" | a hard limit on max-intermediate-size |
All six values are recognized; the path that is optimal in one cost is generally suboptimal in another.
In[169]:=
| OptimalContractionPath |
Out[169]=
{{1,2},{1,2},{2,3},{1,2}}
In[170]:=
| OptimalContractionPath |
Out[170]=
{{1,2},{1,2},{2,3},{1,2}}
In[171]:=
| OptimalContractionPath |
Out[171]=
{{1,2},{3,4},{2,3},{1,2}}
In[172]:=
| OptimalContractionPath |
Out[172]=
{{1,2},{3,4},{2,3},{1,2}}
In[173]:=
| OptimalContractionPath |
Out[173]=
{{1,2},{3,4},{2,3},{1,2}}
"AllowOuterProducts" permits the optimizer to introduce outer-product steps (pairs with no shared index). The default disallows them since they multiply tensor sizes without contracting; some networks benefit from allowing the search to consider them.
In[174]:=
| OptimalContractionPath |
Out[174]=
{{1,2},{1,2},{2,3},{1,2}}
"PreSimplify" mirrors the greedy option – fuses duplicate indices and trivial dimensions before the dynamic-programming search runs.
In[175]:=
| OptimalContractionPath |
Out[175]=
{{1,2},{1,2},{2,3},{1,2}}
4. Path Representations
Three equivalent path representations circulate inside the paclet, with utilities to convert between them. The flat linear form is the opt_einsum list of pairs already seen above; positions are renumbered after every step so each pair refers to the working list at that moment. The tree form is a nested list where each inner list pairs (or groups) two subtrees, recording the contraction order as a binary tree of leaves at the bottom and the final result at the root – convenient for visualization, for engines that consume sub-products in one step, and for any operation that needs to look at sub-contractions as units rather than as positional pairs. and convert in either direction; they are exact inverses up to the choice of how leaves are labeled (numeric tensor positions by default).
In[177]:=
tree=
[tn,Method"flops"]
| PathToTreePath |
| OptimalContractionPath |
Out[177]=
{{3},{{4},{{5},{{1},{2}}}}}
Read the tree from the leaves outward: {1} and {2} pair first (the linear {1,2} step), giving an inner pair {{1},{2}}. Then {5} joins that group, then {4} joins the resulting subtree, and {3} joins last as the outermost pair. The shape is left-leaning because the flops path repeatedly merges the most recently formed intermediate with the next dangling leaf, rather than building two balanced halves; the planner found that this scheduling minimizes total FLOPs given the bond dimensions of this network. TreePathToPath inverts the conversion in the obvious way, recursively flattening the tree into the linear positional form. normalizes a path that may include unary identity steps {i} (which arise when an earlier BinaryTensorNetwork pass binarizes a single hyperedge) by pushing them into pair form. A path is canonical when every step is a pair; both planners return canonical paths.
shows, for each step, which indices were summed out. Pass the path and the hyperedge list. When the network has hyper-edges (any index shared by three or more tensors), use the hyperedges so the index labels match what the planner sees:
In[178]:=
| PathIndexContractions |
| GreedyContractionPath |
| BinaryTensorNetwork |
Out[178]=
{{4},{13,15},{12,14},{7,11}}
Step one sums index 4 (a single shared bond between the first two merged tensors); steps two through four each sum a pair of bonds, including the chain of shared indices that link the two halves of the network after the first merge has happened. The free indices {8, 16} of the network are absent from the list, as expected – they never get contracted, only carried through the intermediates to the final tensor's free legs. This view of the path is the natural input to cost analysis: each summed index contributes a factor of its bond dimension to the FLOP count of the step in which it disappears, so the index-contraction list – together with the size of each index – gives an immediate formula for the path's total cost without re-running the planner.
Just as checks the linear form and checks its stricter pair-only variant, validates the nested-list tree form before it is passed to a contraction engine. A tree is valid when every leaf is a singleton list {i} naming a tensor position and every internal node is a list of two or more subtrees:
In[179]:=
| TreePathQ |
Out[179]=
True
In[180]:=
| TreePathQ |
Out[180]=
True
In[181]:=
| TreePathQ |
5. Visualizing the Contraction
The default labels are symbolic tensor names – each node displays a placeholder T annotated with the dimensions it will carry. This view is useful when you care about which input tensor became which intermediate. To emphasize the dimension flow at every junction, set "Labels" to "Dimensions" instead. Each node then shows just its shape and the contracted axis count, making it easy to spot where intermediates balloon. The ratio between the largest dimension at any node and the dimensions of the leaves is a visual proxy for the path's peak-memory cost:
The tree mirrors what PathToTreePath gave numerically: every internal node is a binary merge of two subtrees, the leaves are the input tensors, and the root is the final result. Reading from leaves to root recovers the exact step-by-step plan that the engine will execute when the expression is activated. The mathematical relationship between the path, the tree, and the final Inactive expression is: the path is a positional schedule, the tree records which tensors meet at which nodes, and the Inactive expression is the same tree decorated with the specific back-end primitives – Inactive[ArrayDot], Inactive[Transpose], and friends – that the chosen Method will execute. Three views, one plan.
6. The Five Contraction Engines
The distinction matters because the same pairwise contraction can be coded several ways, and the cost of the inner kernel – not just the contraction order – affects wall-clock time. Two networks with identical optimal paths will run at very different speeds if one is contracted with "ArrayDot" and the other with "TableSum". A useful mental model: the planner decides what gets contracted together, the engine decides how. To compare two engines, build a small numeric network and contract it twice:
Reproducibility – set a fixed seed before drawing random tensors.
Build a small numeric MPS to benchmark the engines against.
"ArrayDotTranspose" permutes axes to a canonical layout before each pairwise contraction. Same numerical result as "ArrayDot", different internal representation.
Verify against the "ArrayDot" baseline.
"TensorContract" keeps every step as a symbolic TensorContract[TensorProduct[…], …] expression activated at the end. Closest to the textbook definition; slower for large networks.
Same numerics.
Same numerics – the engine is correct, only the end-to-end TensorNetworkContract path is broken.
7. Inactive Expressions
Putting the whole pipeline together: choose a planner, pass its output to TensorNetworkContraction to get an Inactive plan, inspect it (with ContractionTree, with PathIndexContractions, or by reading the symbolic expression directly), then call ActivateTensors when you want a number. For one-shot use, TensorNetworkContract collapses all three steps into a single call. The path is data, the plan is data, the result is data – every level is inspectable and every level is interchangeable, which is why the same paclet handles a three-tensor MPS contraction and a multi-thousand-tensor PEPS overlap with the same vocabulary.
In day-to-day use the most common pattern is the shortest one: TensorNetworkContract[net] by itself. The function selects the greedy heuristic for large networks and the size-optimal dynamic-programming search for small ones, picks "ArrayDot" as the engine, and returns the final tensor. The longer pipeline of this tutorial – choosing the planner, choosing the objective, choosing the engine, holding the plan inactive, walking the tree – is what you reach for when the default takes too long, when memory runs out, when the result must be verified across engines, or when the contraction is part of a symbolic codegen pass that will run thousands of times against tensors of the same shape. The investment pays off precisely when a contraction stops being a one-line call and becomes something you need to reason about.
The "Inactive" option toggles whether the result is an Inactive[…] expression tree (default True, the symbolic form) or the evaluated tensor.
The path argument can also be a method-name string. TensorNetworkContract[tn, "Greedy"] is shorthand for TensorNetworkContract[tn, GreedyContractionPath[tn]]; "Optimal", "flops", "max", "size", "write", "combo", "limit" all dispatch to the corresponding optimal-path call.