Contraction Examples¶

Tensor contraction is the fundamental operation for connecting tensors in a network. It generalizes matrix multiplication to higher-order tensors while automatically enforcing charge conservation at every contracted bond.

Key concepts:

Automatic contraction: Nicole matches indices by tags and opposite directions
Manual specification: Control exactly which indices to contract via position or exclusion
Charge conservation: Only blocks with matching charges on contracted indices contribute
Trace: Special case where a tensor contracts with itself

Contractions in Nicole are symmetry-aware: they only compute blocks that satisfy charge conservation, making them significantly more efficient than dense tensor operations for systems with symmetry.

Basic Contraction¶

group = U1Group()
idx_out = Index(Direction.OUT, group, sectors=(Sector(0, 2), Sector(1, 1)))
idx_in = Index(Direction.IN, group, sectors=(Sector(0, 2), Sector(1, 1)))

# Create two tensors
A = Tensor.random([idx_out, idx_out], itags=["i", "mid"], seed=10)
B = Tensor.random([idx_in, idx_out], itags=["mid", "j"], seed=11)

# Automatic contraction (matching tags with opposite directions)
C = contract(A, B)
print(C)

  info:  2x { 1 x 1 }  having 'A',   Tensor,  { i*, j* }
  data:  2-D float64 (32 B)    2 x 2 => 2 x 2  @ norm = 1.23187

     1.  2x2     |  1x1     [ 0 ; 0 ]    32 B

Manual Pair Specification¶

# Specify axes explicitly: (axis_in_A, axis_in_B)
Cp = contract(A, B, axes=(1, 0))
print(Cp)
print()

# Verify they are identical
diff = (C - Cp).norm()
print(f"Difference from automatic: {diff:.2e}")

  info:  2x { 1 x 1 }  having 'A',   Tensor,  { i*, j* }
  data:  2-D float64 (32 B)    2 x 2 => 2 x 2  @ norm = 1.23187

     1.  2x2     |  1x1     [ 0 ; 0 ]    32 B      

Difference from automatic: 0.00e+00

Matrix-Matrix Multiplication¶

# Two matrices to contract
M1 = Tensor.random([idx_out, idx_in], itags=["i", "j"], seed=20)
M2 = Tensor.random([idx_out, idx_in], itags=["j", "k"], seed=21)

# Multiply: M1 @ M2 (contracts on "j")
result_mm = contract(M1, M2)
print(f"Result indices: {result_mm.itags}")

Result indices: ('i', 'k')

Full Contraction (Scalar Result)¶

# Two tensors that fully contract (need opposite directions for each tag)
T1 = Tensor.random([idx_out, idx_in], itags=["i", "j"], seed=1)
T2 = Tensor.random([idx_in, idx_out], itags=["i", "j"], seed=2)

# Full contraction: all indices contract
scalar = contract(T1, T2)
print(f"Scalar result: {scalar.norm():.4f}")
print(f"Number of indices: {len(scalar.indices)}")

Scalar result: 1.0406
Number of indices: 0

MPS-like Contraction¶

# 3-index tensor: (auxiliary_left, auxiliary_right, physical)
from nicole import conj
idx_aux = Index(Direction.OUT, group, sectors=(Sector(0, 3), Sector(1, 2)))
idx_phys = Index(Direction.OUT, group, sectors=(Sector(0, 1), Sector(1, 1)))

M = Tensor.random([idx_aux, idx_aux.flip(), idx_phys], itags=["left", "right", "phys"], seed=5)
print(f"M tensor:\n{M}\n")

# Contract M with its conjugate on left and phys, excluding right
M_dag = conj(M)
print(f"M† tensor:\n{M_dag}\n")

# Method 1: Specify which axes to contract (axes 0 and 2)
result1 = contract(M, M_dag, axes=((0, 2), (0, 2)))
print(f"Method 1 result (axes=[0,2]):\n{result1}\n")

# Method 2: Exclude right auxiliary index (position 1)
result2 = contract(M, M_dag, excl=((1,),()))
print(f"Method 2 result (excl=1):\n{result2}")

M tensor:

  info:  3x { 3 x 1 }  having 'A',   Tensor,  { left*, right, phys* }
  data:  3-D float64 (152 B)    5 x 5 x 2 => 5 x 5 x 2  @ norm = 3.45509

     1.  3x3x1   |  1x1x1   [ 0 ; 0 ; 0 ]    72 B      
     2.  3x2x1   |  1x1x1   [ 0 ; 1 ; 1 ]    48 B      
     3.  2x2x1   |  1x1x1   [ 1 ; 1 ; 0 ]    32 B      

M† tensor:

  info:  3x { 3 x 1 }  having 'A',   Tensor,  { left, right*, phys }
  data:  3-D float64 (152 B)    5 x 5 x 2 => 5 x 5 x 2  @ norm = 3.45509

     1.  3x3x1   |  1x1x1   [ 0 ; 0 ; 0 ]    72 B      
     2.  3x2x1   |  1x1x1   [ 0 ; 1 ; 1 ]    48 B      
     3.  2x2x1   |  1x1x1   [ 1 ; 1 ; 0 ]    32 B      

Method 1 result (axes=[0,2]):

  info:  2x { 2 x 1 }  having 'A',   Tensor,  { right, right* }
  data:  2-D float64 (104 B)    5 x 5 => 5 x 5  @ norm = 6.40447

     1.  3x3     |  1x1     [ 0 ; 0 ]    72 B      
     2.  2x2     |  1x1     [ 1 ; 1 ]    32 B      

Method 2 result (excl=1):

  info:  2x { 2 x 1 }  having 'A',   Tensor,  { right, right* }
  data:  2-D float64 (104 B)    5 x 5 => 5 x 5  @ norm = 6.40447

     1.  3x3     |  1x1     [ 0 ; 0 ]    72 B      
     2.  2x2     |  1x1     [ 1 ; 1 ]    32 B

Trace¶

# Square tensor
T_trace = Tensor.random([idx_out, idx_in], itags=["i", "i"], seed=99)

# Automatic trace (finds matching itags with opposite directions)
tr = trace(T_trace)
print(f"Trace (automatic): {tr.norm():.4f}")
print(f"Result is scalar: {len(tr.indices) == 0}\n")

# Manual specification by position
tr2 = trace(T_trace, axes=(0, 1))
print(f"Trace (manual): {tr2.norm():.4f}")

Trace (automatic): 0.6905
Result is scalar: True

Trace (manual): 0.6905

Trace Operations¶

# 4-index tensor with paired itags
idx = Index(Direction.OUT, group, sectors=(Sector(0, 2), Sector(1, 1)))
T_multi = Tensor.random(
    [idx, idx.flip(), idx, idx.flip()],
    itags=["a", "a", "b", "b"],
    seed=77
)

# Automatic mode: trace all matching pairs
result_auto = trace(T_multi)
print(f"Result is scalar: {result_auto.is_scalar()}")

# Manual mode: trace specific pair only
partial = trace(T_multi, axes=(0, 1))
print(f"Remaining indices: {partial.itags}")

# Exclusion mode: trace all except specified
partial2 = trace(T_multi, excl=[0, 1])
print(f"Remaining indices: {partial2.itags}")

Result is scalar: True
Remaining indices: ('b', 'b')
Remaining indices: ('a', 'a')

Multiple Contractions¶

# Contract three tensors: A-B-C
A_multi = Tensor.random([idx_out, idx_out], itags=["i", "mid1"], seed=1)
B_multi = Tensor.random([idx_in, idx_out], itags=["mid1", "mid2"], seed=2)
C_multi = Tensor.random([idx_in, idx_out], itags=["mid2", "j"], seed=3)

# Contract step by step
AB = contract(A_multi, B_multi)
print(f"After A-B: {AB.itags}")

ABC = contract(AB, C_multi)
print(f"After A-B-C: {ABC.itags}")

After A-B: ('i', 'mid2')
After A-B-C: ('i', 'j')

Einstein Summation¶

einsum provides a compact subscript-string interface that dispatches to contract, trace, and permute automatically. The -> separator is required; the left-hand side is a comma-separated list of subscripts (one per tensor) and the right-hand side is the output subscript.

Note: Contractions are performed strictly left to right with no optimisation of the contraction order.

Permutation — reorder axes of a single tensor:

T = Tensor.random([idx_out, idx_in], itags=["i", "j"], seed=30)
T_T = einsum('ij->ji', T)
print(f"Original indices: {T.itags}")
print(f"Transposed indices: {T_T.itags}")

Original indices: ('i', 'j')
Transposed indices: ('j', 'i')

Trace to scalar — contract a repeated subscript letter within one tensor:

T_sq = Tensor.random([idx_out, idx_in], itags=["i", "i"], seed=31)
scalar = einsum('ii->', T_sq)
print(f"Trace result is scalar: {scalar.is_scalar()}")
print(f"Value: {scalar.norm():.4f}")

Trace result is scalar: True
Value: 0.3850

Matrix multiply — contract the shared index between two tensors:

A_e = Tensor.random([idx_out, idx_out], itags=["i", "j"], seed=32)
B_e = Tensor.random([idx_in, idx_out], itags=["j", "k"], seed=33)
C_e = einsum('ij,jk->ik', A_e, B_e)
print(f"Result indices: {C_e.itags}")

Result indices: ('i', 'k')

Chain contraction — three tensors contracted left to right:

idx_mid = Index(Direction.OUT, group, sectors=(Sector(0, 2), Sector(1, 1)))
X = Tensor.random([idx_out, idx_mid], itags=["i", "j"], seed=40)
Y = Tensor.random([idx_mid.flip(), idx_mid], itags=["j", "k"], seed=41)
Z = Tensor.random([idx_mid.flip(), idx_in], itags=["k", "l"], seed=42)
result = einsum('ij,jk,kl->il', X, Y, Z)
print(f"Chain result indices: {result.itags}")

Chain result indices: ('i', 'l')