Skip to content

RFC: Distributed Computing Support for Large-Scale Tensor Network Contraction #17

Description

@GiggleLiu

Summary

This RFC proposes adding distributed computing support to omeinsum-rs for large-scale tensor network contraction across compute clusters. As tensor networks grow in size (e.g., quantum circuit simulation, PEPS/MERA for condensed matter), single-node computation becomes infeasible. Distributed execution is essential for:

  • Quantum circuit simulation: 50+ qubit circuits require distributed memory
  • Tensor network machine learning: Large-scale TNML models
  • Condensed matter physics: PEPS/MERA ground state calculations
  • Combinatorial optimization: Large QUBO problems via tropical tensor networks

Current State

What We Have

  • ✅ CPU backend with BLAS acceleration (faer)
  • ✅ CUDA backend with cuTENSOR integration
  • ✅ Contraction order optimization (greedy, TreeSA via omeco)
  • ✅ Multiple algebras (Standard, MaxPlus, MinPlus, MaxMul)
  • ✅ Automatic differentiation for gradients
  • ✅ Stride-based tensors with zero-copy views

What's Missing

  • ❌ Multi-node execution
  • ❌ Distributed tensor storage
  • ❌ Inter-node communication
  • ❌ Distributed contraction scheduling
  • ❌ Fault tolerance

Technology Analysis

Option 1: MPI via mpi crate (rsmpi)

Aspect Assessment
GitHub Stars 577
Downloads/month 4,585
Last Release Dec 2025 (v0.8.1)
Reliability ⭐⭐⭐⭐ (4/5)

Pros:

  • Industry standard for HPC
  • Works on existing cluster infrastructure (SLURM, PBS)
  • High-performance interconnects (InfiniBand, RDMA)
  • Collective operations (allreduce, scatter, gather)
  • Mature ecosystem

Cons:

  • Requires MPI runtime installation
  • Build complexity (needs libclang)
  • Static process model (no dynamic scaling)
  • Error handling is coarse-grained

Option 2: Actor Framework via ractor

Aspect Assessment
GitHub Stars 1,900
Downloads/month 25,796
Last Release Dec 2025 (v0.15.10)
Reliability ⭐⭐⭐⭐ (4/5)

Pros:

  • Pure Rust, no external dependencies
  • Dynamic scaling and fault tolerance
  • Location-transparent messaging
  • Supervision trees for reliability
  • Cloud-native friendly

Cons:

  • Different paradigm from traditional HPC
  • Higher per-message overhead than MPI
  • No native RDMA support
  • Smaller HPC community

Option 3: gRPC via tonic

Aspect Assessment
GitHub Stars 11,800
Downloads/month 9.3M
Last Release Sep 2025 (v0.14.2)
Reliability ⭐⭐⭐⭐⭐ (5/5)

Pros:

  • Extremely mature and widely used
  • Excellent tooling (protobuf, streaming)
  • Works across languages
  • Good for heterogeneous clusters

Cons:

  • Higher latency than MPI
  • Not optimized for HPC workloads
  • Requires service mesh for complex topologies

Option 4: Hybrid (MPI + Rayon)

Aspect Assessment
Rayon Stars 12,600
Rayon Downloads/month 11.4M
Reliability ⭐⭐⭐⭐⭐ (5/5)

Pros:

  • Best of both worlds
  • MPI for inter-node, Rayon for intra-node
  • Standard HPC pattern
  • Optimal resource utilization

Cons:

  • Two parallel models to reason about
  • Thread affinity management complexity

Deep Dive: MPI for Tensor Network Contraction

MPI Communication Patterns for Tensor Networks

Tensor network contraction has specific communication patterns that map well (or poorly) to MPI primitives:

Point-to-Point vs Collective Operations

Operation MPI Primitive Tensor Network Use Case Efficiency
Tensor slice transfer MPI_Send/Recv Moving tensor shards between nodes ⭐⭐⭐ Medium
Broadcast input tensor MPI_Bcast Replicate small tensors to all nodes ⭐⭐⭐⭐⭐ Excellent
Gather partial results MPI_Gather Collect output shards to root ⭐⭐⭐⭐ Good
Reduce contractions MPI_Allreduce Sum partial GEMM results ⭐⭐⭐⭐⭐ Excellent
Scatter input shards MPI_Scatter Distribute tensor along one dim ⭐⭐⭐⭐ Good
All-to-all redistribution MPI_Alltoall Change tensor distribution ⭐⭐ Poor (high cost)

MPI Topology Awareness

MPI can optimize communication based on physical network topology:

// Create Cartesian topology matching cluster layout
let dims = [num_nodes_x, num_nodes_y];
let cart_comm = comm.create_cart(&dims, &[false, false], true);

// Neighbor communication is optimized
cart_comm.shift(0, 1); // Get neighbors in x-direction

Benefits for Tensor Networks:

  • PEPS contractions naturally map to 2D process grids
  • Nearest-neighbor communication dominates
  • Can exploit locality in contraction order

MPI-Specific Distribution Strategies

Strategy A: Replicated Inputs, Sharded Output

Node 0: A (full), B (full) → C[0:n/4, :]
Node 1: A (full), B (full) → C[n/4:n/2, :]
Node 2: A (full), B (full) → C[n/2:3n/4, :]
Node 3: A (full), B (full) → C[3n/4:n, :]

MPI Pattern:

// Broadcast inputs
comm.broadcast(0, &mut tensor_a);
comm.broadcast(0, &mut tensor_b);

// Local computation (embarrassingly parallel)
let local_c = einsum::<A>(&[&a, &b], ...);

// Gather results (optional, if root needs full output)
comm.gather(0, &local_c, &mut full_c);
Pros Cons
Simple to implement Memory inefficient (inputs replicated)
No communication during compute Doesn't scale for large inputs
Perfect load balance Broadcast cost O(n × log(p))

Best for: Small-to-medium inputs, large outputs


Strategy B: SUMMA (Scalable Universal Matrix Multiply Algorithm)

For large matrix contractions ij,jk→ik:

Process Grid (2×2 example):
┌─────────┬─────────┐
│ P(0,0)  │ P(0,1)  │
│ A[0,:], │ A[0,:], │
│ B[:,0]  │ B[:,1]  │
├─────────┼─────────┤
│ P(1,0)  │ P(1,1)  │
│ A[1,:], │ A[1,:], │
│ B[:,0]  │ B[:,1]  │
└─────────┴─────────┘

Each step:
1. Broadcast A block along row
2. Broadcast B block along column
3. Local GEMM
4. Accumulate to C

MPI Pattern:

for k in 0..num_blocks {
    // Broadcast A[i, k] along process row
    let a_owner = (my_row, k % grid_cols);
    row_comm.broadcast(a_owner.1, &mut a_block);
    
    // Broadcast B[k, j] along process column
    let b_owner = (k % grid_rows, my_col);
    col_comm.broadcast(b_owner.0, &mut b_block);
    
    // Local GEMM: C_local += A_block × B_block
    local_gemm(&a_block, &b_block, &mut c_local);
}
Pros Cons
Memory efficient O(n²/p) per node Complex implementation
Communication efficient O(n²/√p) Requires square-ish process grid
Scales to large matrices Block size tuning needed
Overlaps compute and communication Not optimal for small tensors

Best for: Large matrix multiplications, dense tensor contractions


Strategy C: 1D Distribution with Allreduce

For contractions with reduction ij,ij→ (element-wise + sum):

Node 0: A[0:n/4], B[0:n/4] → partial_sum_0
Node 1: A[n/4:n/2], B[n/4:n/2] → partial_sum_1
...
Allreduce(SUM) → total_sum on all nodes

MPI Pattern:

// Each node has a shard
let local_result = local_einsum(&a_shard, &b_shard);

// Reduce across all nodes
let mut global_result = local_result.clone();
comm.allreduce(&local_result, &mut global_result, MPI_SUM);
Pros Cons
Simple, efficient for reductions Only works for full contractions
Allreduce is highly optimized Result replicated (memory overhead)
O(log p) communication rounds

Best for: Trace operations, full tensor contractions, loss computation


Strategy D: Tree-Based Pairwise Contraction

For multi-tensor networks A,B,C,D,E → result:

Contraction Tree (from omeco optimizer):
        result
       /      \
     AB        CDE
    /  \      /   \
   A    B    CD    E
            /  \
           C    D

Distributed Assignment:
- Nodes 0-1: Compute A⊗B → AB
- Nodes 2-3: Compute C⊗D → CD  
- Node 2: Compute CD⊗E → CDE (after receiving CD)
- Node 0: Compute AB⊗CDE → result (after receiving both)

MPI Pattern:

match my_assignment {
    Task::Contract(a, b, dest) => {
        let result = local_einsum(&tensors[a], &tensors[b]);
        if dest != my_rank {
            comm.send(dest, &result);
        }
    }
    Task::Receive(src, then_contract) => {
        let received = comm.recv(src);
        // Continue with next contraction
    }
}
Pros Cons
Exploits contraction order optimization Load imbalance (tree structure)
Minimizes intermediate tensor size Complex scheduling
Natural for tensor networks Point-to-point latency sensitive

Best for: Complex tensor networks with many tensors


MPI Performance Considerations

1. Latency vs Bandwidth

Interconnect Latency Bandwidth Best Strategy
InfiniBand EDR ~1 μs 100 Gb/s Fine-grained, SUMMA
InfiniBand HDR ~0.6 μs 200 Gb/s Fine-grained, SUMMA
100G Ethernet ~5 μs 100 Gb/s Coarse-grained, fewer messages
10G Ethernet ~50 μs 10 Gb/s Very coarse, batch transfers

Rule of thumb: Message size > bandwidth × latency for efficiency

  • InfiniBand: Messages > 100 KB are efficient
  • Ethernet: Messages > 500 KB preferred

2. MPI Implementation Differences

Implementation Strengths Considerations
OpenMPI Flexible, many transports Default for most clusters
MPICH Reference impl, ABI standard Good compatibility
Intel MPI Optimized for Intel hardware Best on Intel clusters
MVAPICH2 Best InfiniBand support HPC-focused
Cray MPI Cray system optimized Supercomputers only

The mpi crate (rsmpi) works with all of these.

3. CUDA-Aware MPI

For GPU clusters, CUDA-aware MPI eliminates host staging:

// Without CUDA-aware MPI:
gpu_tensor.copy_to_host(&mut host_buffer);
comm.send(dest, &host_buffer);  // CPU memory
// ... on receiver ...
comm.recv(src, &mut host_buffer);
gpu_tensor.copy_from_host(&host_buffer);

// With CUDA-aware MPI:
comm.send(dest, gpu_tensor.as_ptr());  // Direct GPU memory!
comm.recv(src, gpu_tensor.as_mut_ptr());

Supported by: OpenMPI, MVAPICH2, Cray MPI
Speedup: 2-10x for GPU-heavy workloads

4. Non-Blocking Operations

Overlap communication with computation:

// Start async receive for next iteration
let recv_req = comm.irecv(src, &mut next_buffer);

// Compute with current data
let result = local_einsum(&current_a, &current_b);

// Wait for communication to complete
recv_req.wait();

MPI Challenges for Tensor Networks

Challenge 1: Dynamic Tensor Shapes

MPI typically requires knowing message sizes upfront:

// Problem: receiver doesn't know tensor shape
// Solution 1: Send shape first
comm.send(dest, &tensor.shape());  // Small metadata message
comm.send(dest, tensor.data());    // Large data message

// Solution 2: Use MPI_Probe
let status = comm.probe(src);
let count = status.count();  // Learn size before receiving

Challenge 2: Irregular Communication

Tensor networks often have irregular patterns unlike regular stencils:

// PEPS boundary contraction - neighbors only
for neighbor in tensor_network.neighbors(my_region) {
    comm.sendrecv(neighbor, &boundary_data, &mut incoming);
}

Mitigation: Use MPI_Neighbor_alltoall with graph topology

Challenge 3: Load Balancing

Contraction costs vary dramatically:

// Tensor dimensions affect cost cubically
let cost_a = a.shape().iter().product::<usize>();  // e.g., 1000
let cost_b = b.shape().iter().product::<usize>();  // e.g., 1000000

// Static assignment leads to imbalance
// Solution: Work stealing or dynamic scheduling

Challenge 4: Memory Constraints

Large intermediate tensors may exceed node memory:

// Check before allocation
let required = intermediate_size * size_of::<f64>();
if required > available_memory {
    // Use tiled/streaming algorithm
    // Or spill to disk/NVMe
}

Comparison: MPI vs Alternatives for Tensor Networks

Aspect MPI Actors (ractor) gRPC (tonic)
Latency ~1 μs (IB) ~100 μs ~1 ms
Bandwidth 200 Gb/s Network limited Network limited
RDMA ✅ Native
GPU-Direct ✅ CUDA-aware
Fault Tolerance ❌ Process fails = job fails ✅ Supervision ✅ Retries
Dynamic Scaling ❌ Static ✅ Runtime spawn ✅ Service mesh
Debugging ⚠️ Complex ✅ Better tooling ✅ Standard tooling
Cloud Deploy ⚠️ Needs MPI runtime ✅ Native ✅ Native
HPC Cluster ✅ Standard ⚠️ Unusual ⚠️ Unusual

Recommended Approach: Hybrid MPI + Rayon

Based on analysis, we recommend a hybrid approach with abstraction layers:

┌─────────────────────────────────────────────────────────────┐
│                     User API Layer                          │
│  einsum_distributed<A, T>(&tensors, &indices, &output)     │
└─────────────────────────────────────────────────────────────┘
                              │
                              ▼
┌─────────────────────────────────────────────────────────────┐
│              Distributed Scheduler                          │
│  - Partition contraction tree across nodes                  │
│  - Schedule data movement                                   │
│  - Handle load balancing                                    │
└─────────────────────────────────────────────────────────────┘
                              │
              ┌───────────────┼───────────────┐
              ▼               ▼               ▼
┌─────────────────┐ ┌─────────────────┐ ┌─────────────────┐
│   Node 0        │ │   Node 1        │ │   Node N        │
│ ┌─────────────┐ │ │ ┌─────────────┐ │ │ ┌─────────────┐ │
│ │Local Einsum │ │ │ │Local Einsum │ │ │ │Local Einsum │ │
│ │ (CPU/GPU)   │ │ │ │ (CPU/GPU)   │ │ │ │ (CPU/GPU)   │ │
│ └─────────────┘ │ │ └─────────────┘ │ │ └─────────────┘ │
│ ┌─────────────┐ │ │ ┌─────────────┐ │ │ ┌─────────────┐ │
│ │   Rayon     │ │ │ │   Rayon     │ │ │ │   Rayon     │ │
│ │ (intra-node)│ │ │ │ (intra-node)│ │ │ │ (intra-node)│ │
│ └─────────────┘ │ │ └─────────────┘ │ │ └─────────────┘ │
└────────┬────────┘ └────────┬────────┘ └────────┬────────┘
         │                   │                   │
         └───────────────────┼───────────────────┘
                             │
                    ┌────────▼────────┐
                    │  Communication  │
                    │     Layer       │
                    │  (MPI / Actor)  │
                    └─────────────────┘

Tiered Backend Selection

pub enum CommunicationBackend {
    /// MPI for HPC clusters (best performance)
    Mpi(MpiCommunicator),
    
    /// Actors for cloud/fault-tolerant (good flexibility)
    Actor(RactorCommunicator),
    
    /// gRPC for heterogeneous systems (best compatibility)
    Grpc(TonicCommunicator),
    
    /// Single-node with Rayon (baseline)
    Local(RayonExecutor),
}

impl CommunicationBackend {
    /// Auto-detect best backend for environment
    pub fn detect() -> Self {
        if mpi_available() && on_hpc_cluster() {
            Self::Mpi(MpiCommunicator::new())
        } else if kubernetes_detected() {
            Self::Actor(RactorCommunicator::new())
        } else {
            Self::Local(RayonExecutor::new())
        }
    }
}

Priority Order for Implementation:

  1. MPI - Primary target for HPC users (largest performance gains)
  2. Rayon - Intra-node parallelism (complements MPI)
  3. Actors - Cloud deployment and fault tolerance
  4. gRPC - Heterogeneous and polyglot environments

Proposed Architecture

1. Communication Trait Abstraction

/// Abstract communication layer - can be MPI, actors, or gRPC
pub trait Communicator: Clone + Send + Sync {
    type Error: std::error::Error;
    
    fn rank(&self) -> usize;
    fn size(&self) -> usize;
    
    // Point-to-point
    async fn send<T: Scalar>(&self, dest: usize, data: &[T]) -> Result<(), Self::Error>;
    async fn recv<T: Scalar>(&self, src: usize, len: usize) -> Result<Vec<T>, Self::Error>;
    
    // Collective operations
    async fn broadcast<T: Scalar>(&self, root: usize, data: &mut [T]) -> Result<(), Self::Error>;
    async fn allreduce<T: Scalar, F: Fn(T, T) -> T>(&self, data: &mut [T], op: F) -> Result<(), Self::Error>;
    async fn scatter<T: Scalar>(&self, root: usize, send: &[T], recv: &mut [T]) -> Result<(), Self::Error>;
    async fn gather<T: Scalar>(&self, root: usize, send: &[T], recv: &mut [T]) -> Result<(), Self::Error>;
    
    // Synchronization
    async fn barrier(&self) -> Result<(), Self::Error>;
}

2. Distributed Tensor

/// Tensor distributed across multiple nodes
pub struct DistributedTensor<T: Scalar, B: Backend, C: Communicator> {
    /// Local shard of the tensor
    local: Tensor<T, B>,
    /// Distribution strategy
    distribution: Distribution,
    /// Global shape
    global_shape: Vec<usize>,
    /// Communicator handle
    comm: C,
}

pub enum Distribution {
    /// Tensor replicated on all nodes
    Replicated,
    /// Tensor sharded along specific dimension
    Sharded { dim: usize, chunks: Vec<(usize, usize)> },
    /// Block-cyclic distribution (ScaLAPACK style)
    BlockCyclic { block_size: Vec<usize> },
}

3. Distributed Contraction Strategies

pub enum ContractionStrategy {
    /// Each node computes part of output (embarrassingly parallel)
    OutputParallel {
        output_partition: Vec<Range<usize>>,
    },
    
    /// Tensor A sharded, B replicated, reduce partial results
    ReduceScatter {
        shard_dim: usize,
    },
    
    /// SUMMA-style 2D distribution for large matrix multiply
    Summa {
        proc_grid: (usize, usize),
        block_size: (usize, usize),
    },
    
    /// Tree-based reduction for multi-tensor contractions
    TreeReduction {
        tree: NestedEinsum<usize>,
        node_assignment: Vec<usize>,
    },
}

4. Distributed Einsum API

impl<L: Label> Einsum<L> {
    /// Optimize for distributed execution
    pub fn optimize_distributed(&mut self, comm: &impl Communicator) -> &mut Self {
        // Consider communication costs in optimization
        // Assign subtrees to nodes
        // Minimize data movement
    }
    
    /// Execute across cluster
    pub async fn execute_distributed<A, T, B, C>(
        &self,
        tensors: &[DistributedTensor<T, B, C>],
        comm: &C,
    ) -> Result<DistributedTensor<T, B, C>, DistributedError>
    where
        A: Algebra<Scalar = T>,
        T: Scalar,
        B: Backend,
        C: Communicator,
    {
        // 1. Analyze data locality
        // 2. Schedule contractions to minimize communication
        // 3. Execute local contractions in parallel
        // 4. Communicate intermediate results
        // 5. Combine final result
    }
}

Implementation Phases

Phase 1: Foundation (Core Infrastructure)

  • Define Communicator trait
  • Implement MPI backend (MpiCommunicator)
  • Implement mock communicator for testing
  • Add DistributedTensor type
  • Basic point-to-point tensor transfer
  • Serialization for tensor data

Phase 2: Basic Distribution

  • Implement Distribution strategies
  • Tensor sharding and gathering
  • Replicated tensor broadcast
  • Simple output-parallel contraction
  • Integration tests on multi-process

Phase 3: Optimized Contractions

  • Communication-aware contraction optimizer
  • SUMMA-style distributed GEMM
  • Tree-based reduction for multi-tensor
  • Overlap computation and communication
  • Load balancing heuristics

Phase 4: Advanced Features

  • Rayon integration for intra-node parallelism
  • CUDA-aware MPI for GPU clusters
  • Checkpointing for long-running contractions
  • Fault tolerance (actor-based fallback)
  • Dynamic load balancing

Phase 5: Ecosystem Integration

  • Actor-based communicator (RactorCommunicator)
  • Cloud deployment support (Kubernetes)
  • Benchmark suite for distributed performance
  • Documentation and examples

Challenges and Considerations

1. Communication Cost Modeling

The contraction optimizer (omeco) currently optimizes for FLOPs. For distributed execution, we need:

/// Cost model for distributed contraction
pub struct DistributedCost {
    /// Compute cost (FLOPs)
    compute: f64,
    /// Communication cost (bytes transferred)
    communication: f64,
    /// Memory high-water mark per node
    memory_per_node: f64,
    /// Load imbalance factor
    imbalance: f64,
}

2. Tensor Redistribution

When contraction order requires different distributions:

// A is sharded on dim 0, B is sharded on dim 1
// For A @ B, we need to redistribute
async fn redistribute<T, B, C>(
    tensor: &DistributedTensor<T, B, C>,
    new_dist: Distribution,
) -> DistributedTensor<T, B, C>

3. Memory Management

Large tensors may not fit in single-node memory:

  • Streaming/tiled execution
  • Out-of-core algorithms
  • Memory-aware scheduling

4. Heterogeneous Clusters

Mixed CPU/GPU nodes:

pub struct NodeCapabilities {
    cpus: usize,
    gpus: Vec<GpuInfo>,
    memory: usize,
    interconnect: Interconnect,
}

5. Debugging and Profiling

Distributed debugging is hard:

  • Per-node logging with timestamps
  • Communication tracing
  • Deadlock detection
  • Performance profiling (communication vs compute)

References

  1. Cyclops Tensor Framework - Distributed tensor contractions in C++
  2. ITensor - Tensor network library with MPI support
  3. ExaTENSOR - Exascale tensor algebra
  4. cuQuantum - NVIDIA's distributed quantum simulation
  5. opt_einsum - Python einsum optimizer (for algorithm reference)

Questions for Discussion

  1. MPI vs Actors: Should we prioritize MPI (HPC clusters) or actors (cloud-native) first?
  2. API Design: Should distributed execution be opt-in or should the API automatically detect cluster availability?
  3. Scope: Should we support arbitrary tensor distributions or focus on specific patterns (sharded, replicated)?
  4. Integration: How should this interact with the existing Backend trait?

/cc @maintainers

🤖 Generated with Claude Code

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions