Einsum-Based Multiplicative Update Framework

• The einsum-based multiplicative update framework is a unifying approach for nonnegative tensor factorization using efficient tensor contraction and majorization–minimization.
• It leverages flexible model strings to specify arbitrary contraction patterns and accommodates a broad spectrum of loss functions including the $(α,β)$-divergence family.
• Empirical results show significant acceleration and improved interpretability over traditional methods, especially in handling large, sparse, or incomplete tensor data.

The einsum-based multiplicative update framework is a unifying approach for fitting a broad class of nonnegative tensor factorizations using majorization–minimization (MM) and accelerated by efficient tensor contraction (“einsum”) operations. This framework generalizes standard CP and Tucker multiplicative updates to any model expressible as a tensor contraction, accommodates a spectrum of loss functions including the full $(\alpha, \beta)$-divergence family, and scales to extremely large, sparse, or incomplete tensors. It provides a modeling interface that allows users to specify arbitrary contraction patterns via simple model strings, facilitating flexible, domain-tailored decompositions without manual derivation of update rules (Hood et al., 2 Feb 2026, Leplat, 16 Feb 2026).

1. Mathematical Foundations and Model Generality

Let $$Y \in \mathbb{R}_{\ge0}^{I_1 \times \cdots \times I_M}$$ denote an $M$-mode observed nonnegative tensor. The objective is to find a nonnegative approximation $\widehat Y$ minimizing an entrywise sum of divergence losses: $$L(Y, \widehat Y) = \sum_{i_1,\dots,i_M} d(y_{i_1\cdots i_M}, \widehat y_{i_1\cdots i_M}),$$ subject to $\widehat y_{i_1\cdots i_M} \ge 0$.

The factorization is described by a set of nonnegative tensors $\Theta^{(1)},\ldots,\Theta^{(L)}$, each affiliated with observed and latent index sets. The generative model for tensor entries is formulated as

$$\widehat y_{\text{obs}} = \sum_{\text{lat}} \prod_{\ell=1}^L \theta^{(\ell)}_{\text{obs}_\ell, \text{lat}_\ell},$$

where $\text{obs}$ are observed indices, $\text{lat}$ are latent (summed) indices, and $$Y \in \mathbb{R}_{\ge0}^{I_1 \times \cdots \times I_M}$$0 are the index sets for $$Y \in \mathbb{R}_{\ge0}^{I_1 \times \cdots \times I_M}$$1. This structure subsumes standard decompositions (CP, Tucker, tensor-train) and custom contraction models.

The einsum framework specifies this via a string, e.g., "i r, j r, k r -> i j k" for a 3-way CP factorization, allowing automated construction of contraction and update routines (Hood et al., 2 Feb 2026).

2. Loss Functions and $$Y \in \mathbb{R}_{\ge0}^{I_1 \times \cdots \times I_M}$$2-Divergences

A key feature is support for a wide class of loss functions: $$Y \in \mathbb{R}_{\ge0}^{I_1 \times \cdots \times I_M}$$3 with the $$Y \in \mathbb{R}_{\ge0}^{I_1 \times \cdots \times I_M}$$4-divergence defined as: $$Y \in \mathbb{R}_{\ge0}^{I_1 \times \cdots \times I_M}$$5 This formulation unifies Frobenius (squared-Euclid), Kullback–Leibler, Itakura–Saito, and other divergences, enabling choice according to the generative assumptions or application (Hood et al., 2 Feb 2026, Leplat, 16 Feb 2026).

3. Multiplicative Update Derivation via Majorization–Minimization

The MM methodology constructs a tight surrogate for the loss over current factors, then analytically minimizes it for each block update:

• The divergence is decomposed into convex and concave parts with respect to $$Y \in \mathbb{R}_{\ge0}^{I_1 \times \cdots \times I_M}$$6.
• Using the essential “decomposability property,” the resulting update for each $$Y \in \mathbb{R}_{\ge0}^{I_1 \times \cdots \times I_M}$$7 takes the multiplicative form: $$Y \in \mathbb{R}_{\ge0}^{I_1 \times \cdots \times I_M}$$8 where $$Y \in \mathbb{R}_{\ge0}^{I_1 \times \cdots \times I_M}$$9 are contractions (einsum over factors and auxiliary tensors $M$0, $M$1 derived from the loss), $M$2 is invertible and loss-dependent, and einsum strings are constructed by swapping the current factor’s indices with model outputs (Hood et al., 2 Feb 2026).

For CP and Tucker models under $M$3-divergence, these updates reduce to repeated einsum operations with powers and masks, as formalized in the block-MM and joint-MM pseudocode (Leplat, 16 Feb 2026).

Model	Update Formulation	Einsum Example (CP, 3-way)
Block-MM	$M$4	einsum('ijk,jr,kr->ir', P, A2, A3)
Joint-MM	Reference-cached surrogates, inner updates	Uses $M$5

4. Algorithmic Implementation and Computational Aspects

The implementation involves the following steps:

1. User specifies $M$6, a model string (model_str), $M$7 (divergence), and initial $M$8.
2. Precompute einsum strings for each factor.
3. Iterate until convergence: a. Compute reconstruction $M$9. b. For each factor, compute numerators and denominators using auxiliary functions $\widehat Y$0 and einsum contractions. c. Apply the block update formula.
4. On completion, factors are returned; loss, iterations, and wall time can be monitored (Hood et al., 2 Feb 2026).

Efficient implementation leverages hardware-optimized einsum routines (e.g., via BLAS or GPU backends), yielding per-iteration costs linear in the number of nonzeros plus model parameters. In practice, the method scales to tensors with hundreds of millions of entries in seconds.

Empirical complexity benchmarks demonstrate that einsum-based updates typically outperform classical unfolding-based routines (MTTKRP/Khatri-Rao), with speedups up to $\widehat Y$1 in single-thread benchmarks and even greater acceleration when parallelized (Leplat, 16 Feb 2026).

5. Treatment of Missing Data

Handling of missing data is implemented by introducing a binary mask $\widehat Y$2: $\widehat Y$3 in all loss and update computations. Masked-out entries are omitted from all surrogate contractions, efficiently excluding them from the optimization (Hood et al., 2 Feb 2026).

This mechanism easily integrates with the einsum-based contraction pattern, requiring only one additional masking operation with negligible computational overhead.

6. Interface, Model Specification, and Extendability

The framework provides a user interface where models are specified by einsum-compatible strings. The user selects or defines the loss (e.g., via $\widehat Y$4), and optionally custom-initializes factor tensors. Batch runs, convergence settings, restarts, and validation-based early stopping are supported.

Beyond standard CP and Tucker factorizations, this enables custom or hybrid models—e.g., non-halting, block-structured, domain-aware contractions—without solver modification. Loss functions can be extended, provided that the convex–concave decomposability and invertibility conditions are met. Extensions include:

• Discrete and count losses (Bernoulli, negative binomial)
• Structured sparsity via explicit $\widehat Y$5 or $\widehat Y$6 penalties
• Probabilistic or Bayesian surrogates (Polya-urn MM updates)
• Hybridization with automatic differentiation for differentiable submodules (Hood et al., 2 Feb 2026)

7. Empirical Performance and Limitations

On large real-world tensors such as spatio-temporal Uber pickups, ICEWS event counts, and WITS trade flows, einsum-based MU (notably, the NNEinFact package) consistently converged up to $\widehat Y$7 faster than Adam and standard gradient routines, achieved over $\widehat Y$8 reductions in held-out divergence relative to off-the-shelf CP/Tucker models, and produced more interpretable latent components (e.g., morning commutes, nightlife). Qualitative patterns are more sharply resolved with custom domain-encoded contractions (Hood et al., 2 Feb 2026).

Limitations include:

• Requirement that loss admits a suitable convex–concave, decomposable surrogate.
• Convergence only to stationary points or blockwise local minima, with sensitivity to initialization; multiple restarts are recommended.
• Failure or divergence if auxiliary ratios $\widehat Y$9 escape the invertibility domain of $$L(Y, \widehat Y) = \sum_{i_1,\dots,i_M} d(y_{i_1\cdots i_M}, \widehat y_{i_1\cdots i_M}),$$0.
• For the joint-MM method, while each inner step decreases the surrogate, asymptotic convergence of the joint sequence is not guaranteed in all cases (Leplat, 16 Feb 2026).

A plausible implication is that the einsum-based multiplicative update paradigm will subsume many bespoke methods in nonnegative multiway modeling, provided loss decomposability and contraction structures are carefully verified.

• "Near-Universal Multiplicative Updates for Nonnegative Einsum Factorization" (Hood et al., 2 Feb 2026)
• "Joint Majorization-Minimization for Nonnegative CP and Tucker Decompositions under $$L(Y, \widehat Y) = \sum_{i_1,\dots,i_M} d(y_{i_1\cdots i_M}, \widehat y_{i_1\cdots i_M}),$$1-Divergences: Unfolding-Free Updates" (Leplat, 16 Feb 2026)