---
title: Einsum-Based Multiplicative Update
url: https://www.emergentmind.com/topics/einsum-based-multiplicative-update
type: topic
---

# Einsum-Based Multiplicative Update

An einsum-based multiplicative update is a general framework for fitting nonnegative tensor factorizations by casting the factorization as a sequence of tensor contractions parameterized by Einstein summation (einsum) notation. This approach, as implemented in NNEinFact, enables the application of multiplicative update methods to essentially any nonnegative tensor factorization model that can be expressed as a tensor contraction, with a broad variety of loss functions including the $(α,β)$-divergence. The algorithm utilizes Python-style einsum strings to specify custom factorization models, combines a universal update formula rooted in majorization–minimization theory, and offers practical routines for scaling, handling missing data, and GPU acceleration [2602.02759].

## 1. General Framework: Tensor Factorization as Einsum

Nonnegative tensor factorization aims to approximate a nonnegative $M$-order tensor $Y \in \mathbb{R}_+^{I_1 \times \cdots \times I_M}$ with a parameterized factorization $\hat Y$ constructed via $L$ factor tensors. The approximation is given by an $L$-way tensor contraction:
$$
\hat y_{\mathcal{i}} = \sum_{r_1=1}^{R_1} \cdots \sum_{r_K=1}^{R_K} \prod_{\ell=1}^L \Theta^{(\ell)}_{i_\ell, r_\ell}
$$
where each factor $\Theta^{(\ell)}$ is a tensor whose modes correspond to subsets of the observed and latent indices. The contraction pattern is specified by an einsum string, for example,
```
model_str = "i r1, j r1, a r2, r1 r2 -> i j a"
```
invokes the contraction to yield the desired shape and structure of $\hat Y$. Custom models, including CP, Tucker, tensor train, or arbitrary user-defined forms, can be defined by varying the einsum string [2602.02759].

## 2. Universal Multiplicative Update Derivation

The general goal is to minimize a differentiable loss
$$
L(Y,\hat Y) = \sum_{\text{indices }i} \ell(y_i, \hat y_i)
$$
where each factor $\Theta^{(\ell)} \geq 0$. The core update is a multiplicative step for each factor $\ell$:
$$
\Theta^{(\ell)} \leftarrow \Theta^{(\ell)} \odot g^{-1}\left(\frac{\text{Numerator}}{\text{Denominator}}\right)
$$
where the numerator and denominator are defined via sum contractions over all relevant indices involving element-wise functions $a(\cdot,\cdot)$ and $b(\cdot,\cdot)$, and $g(\cdot)$ is a scalar map dictated by the loss choice. Specifically,
- $\text{Numerator} = \sum_{\text{sum-indices}} a(Y,\hat Y)\times$ (other factors)
- $\text{Denominator} = \sum_{\text{sum-indices}} b(Y,\hat Y)\times$ (other factors)

Alternatively, in the positive/negative gradient view:
$$
\Theta^{(\ell)}_{\text{new}} = \Theta^{(\ell)}_{\text{old}} \odot \frac{[\nabla^-_\ell L]}{[\nabla^+_\ell L]}
$$
with contractions corresponding to appropriate einsum patterns for the positive and negative components [2602.02759].

Common loss functions and their corresponding $a$, $b$, $g$ mappings are:

| Loss type        | $a(y,\hat y)$           | $b(y,\hat y)$         | $g(\lambda)$   |
|------------------|------------------------|-----------------------|---------------|
| Frobenius        | $y$                    | $\hat y$              | $\lambda$     |
| KL Divergence    | $y/\hat y$             | $1$                   | $\lambda$     |
| $(α,β)$-Divergence | $y^α \hat y^{β-1}$   | $\hat y^{α+β-1}$      | $\lambda^{|1-β|}$ (regime-dependent) |

For each loss, the corresponding update embodies the majorization–minimization step derived under mild convexity and decomposability conditions [2602.02759].

## 3. Computational Procedure and Pseudocode

The practical algorithm cycles through the following steps until convergence:

1. **Recompute Fit:** Compute $\hat Y$ via einsum contraction.
2. **Per-Factor Update:** For each factor $\ell$:
   - Compute the numerator (A): einsum contraction of all factors with $a(Y,\hat Y)$ at position $\ell$
   - Compute the denominator (B): einsum contraction with $b(Y,\hat Y)$
   - Update: $\Theta^{(\ell)} \leftarrow \Theta^{(\ell)} \odot g^{-1}(A/B)$
3. **Convergence Check:** Early stopping via held-out divergence (5–10%) or relative loss decrease $< 10^{-6}$

Python-style pseudocode:

```python
# Precompute per‐factor swap of output index
einstr = [swap(model_str, l) for l in range(L)]
while not converged:
    Yhat = einsum(model_str, Θ[0],…,Θ[L−1])
    for l in range(L):
        A = einsum(einstr[l], Θ[0],…,Θ[l−1], a(Y,Yhat), Θ[l+1],…,Θ[L−1])
        B = einsum(einstr[l], Θ[0],…,Θ[l−1], b(Y,Yhat), Θ[l+1],…,Θ[L−1])
        Θ[l] *= g_inv(A / B)
```
This routine allows rapid prototyping of new factorization models by simply specifying a new einsum string and loss function [2602.02759].

## 4. Theoretical Guarantees and Convergence

The update scheme follows a majorization–minimization (MM) framework. At each block update, a tight surrogate $Q(\Theta^{(\ell)}|\Theta^{(\ell)}_{\text{old}}) \geq L$ is constructed, with equality at the current iterate. Minimizing $Q$ with respect to $\Theta^{(\ell)}$ yields the stated multiplicative update. Theorem 2.3 asserts that, assuming decomposability of the loss and model structure, each update strictly decreases the loss $L$ and the sequence converges to a stationary point (local minimum) [2602.02759]. This holds for a wide class of models and losses supported by the einsum formulation.

## 5. Computational Complexity and Implementation Notes

The per-iteration complexity per factor $\ell$ consists of two einsum computations of cost $O(\prod_{\ell' \neq \ell} I_{\ell'} R_{\ell'})$, and an $O(I_\ell R_\ell)$ cost to apply the update. For sparse $Y$, computation can be restricted to nonzero entries, reducing cost to $O(\operatorname{nnz}(Y)\sum R_\ell)$. Only $L$ factor matrices/tensors and the current $\hat Y$ are stored in memory. Missing data can be incorporated by maintaining a mask $M \in \{0,1\}^{I_1\cdots I_M}$ and replacing $Y\to M\odot Y$, $\hat Y\to M\odot\hat Y$ in all contractions. GPU-accelerated einsum (via PyTorch or NumPy backends) yields multifold speedups relative to explicit loops [2602.02759].

## 6. Empirical Performance and Applications

Empirical evaluation considers large, real-world datasets, such as:
- Uber (5-way, $27\times 7\times 24\times 100\times 100$)
- ICEWS (4-way, $249^2\times 20\times 228$)
- WITS (4-way, $196^2\times 96\times 29$)

Optimization comparisons demonstrate that NNEinFact outperforms Adam (tested over six learning rates) across $(α,β)$ losses, converging up to 90$\times$ faster in wall time and achieving less than half the held-out divergence of gradient-based methods. Model-level comparisons show that custom einsum models fitted by this approach yield 10–37% lower held-out divergence than standard CP, Tucker, tensor train, or low-rank Tucker models. Qualitative analyses (e.g., the Uber dataset) recover interpretable spatiotemporal classes using significantly fewer parameters relative to CP models [2602.02759].

## 7. Applicability and Significance

The einsum-based multiplicative update framework provides a universal, efficient, and customizable method for nonnegative tensor factorization across a diversity of losses and model structures. It enables researchers to flexibly define bespoke tensor contractions and losses, eliminates the need to design new update rules for each model, and empirically achieves superior convergence and predictive performance relative to standard approaches. This framework supports missing data, scales to tensors with hundreds of millions of entries, and is amenable to GPU acceleration, facilitating practical deployment across scientific domains with complex multiway data [2602.02759].

Source: https://www.emergentmind.com/topics/einsum-based-multiplicative-update