MTTKRP in CP Tensor Decomposition

• MTTKRP is a key tensor operation that multiplies a tensor unfolding with the Khatri-Rao product of remaining factor matrices.
• It plays a fundamental role in CP decomposition, serving as the gradient-bearing step in alternating least-squares and generalized optimization methods.
• MTTKRP performance is challenged by memory-bound, irregular sparse cases and communication-intensive dense cases, spurring algorithmic and hardware innovations.

Searching arXiv for recent and foundational papers on MTTKRP to ground the article. Searching arXiv for "Matricized Tensor Times Khatri-Rao Product MTTKRP". Matricized Tensor Times Khatri-Rao Product (MTTKRP) is the mode-wise tensor kernel that multiplies a tensor unfolding by the Khatri-Rao product of the remaining factor matrices. It is the canonical workhorse of CANDECOMP/PARAFAC (CP) and CP/PARAFAC tensor decomposition, appears in alternating least squares (ALS), gradient-based methods, generalized CP objectives, and constrained variants, and is repeatedly identified as the dominant computational bottleneck. In sparse regimes, the kernel is characterized as memory-bound, irregular, and synchronization-heavy; in dense regimes, its performance is strongly shaped by communication cost and data layout rather than arithmetic alone (Wijeratne et al., 20 Jul 2025, Ballard et al., 2017).

1. Algebraic definition and equivalent forms

For an $N$-mode tensor $\mathcal{X}$, CP decomposition seeks factor matrices

$F_0, F_1, \ldots, F_{N-1}$

such that

$$\mathcal{X} \approx [\![\lambda; F_0, F_1, \ldots, F_{N-1}]\!] = \sum_{r=0}^{R-1} \lambda_r \, f^{(0)}_r \circ f^{(1)}_r \circ \cdots \circ f^{(N-1)}_r.$$

Within this model, the mode-$n$ MTTKRP is

$$X_{(n)} \left( F_{N-1} \odot \cdots \odot F_{n+1} \odot F_{n-1} \odot \cdots \odot F_0 \right),$$

where $X_{(n)}$ is the mode-$n$ unfolding and $\odot$ is the Khatri-Rao product (Wijeratne et al., 20 Jul 2025).

The Khatri-Rao product is the column-wise Kronecker product. If $A \in \mathbb{R}^{I \times R}$ and $\mathcal{X}$0, then

$\mathcal{X}$1

and its $\mathcal{X}$2-th column is $\mathcal{X}$3. The same construction extends across all input modes except the output mode. This column-wise structure is what couples factor matrices to the tensor data in CP updates (Wijeratne et al., 20 Jul 2025).

In sparse form, MTTKRP is most naturally expressed as an accumulation over nonzeros. For a nonzero tensor element

$\mathcal{X}$4

the contribution to the output factor row for mode $\mathcal{X}$5 is

$\mathcal{X}$6

For a third-order tensor, the same structure appears in the familiar form

$\mathcal{X}$7

or, in COO iteration form,

$\mathcal{X}$8

These equations make explicit that each nonzero induces a rank-$\mathcal{X}$9 vector update into one output row (Wijeratne et al., 2021).

The same kernel admits a cost distinction between dense and sparse derivative or data tensors. In the generalized CP setting, the mode-$F_0, F_1, \ldots, F_{N-1}$0 MTTKRP

$F_0, F_1, \ldots, F_{N-1}$1

has cost $F_0, F_1, \ldots, F_{N-1}$2 for dense $F_0, F_1, \ldots, F_{N-1}$3 and $F_0, F_1, \ldots, F_{N-1}$4 for sparse $F_0, F_1, \ldots, F_{N-1}$5. This dense-versus-sparse dichotomy recurs throughout the literature: the algebraic definition is unchanged, but the dominant implementation issues shift from BLAS organization to irregular access, synchronization, and communication (Kolda et al., 2019).

2. Role in CP, generalized CP, and constrained optimization

In least-squares CP, MTTKRP is not merely an auxiliary tensor primitive; it is the gradient-bearing data term. For the objective

$F_0, F_1, \ldots, F_{N-1}$6

the gradient with respect to factor $F_0, F_1, \ldots, F_{N-1}$7 is

$F_0, F_1, \ldots, F_{N-1}$8

This is why the kernel appears in ALS, Gauss-Newton, damped Gauss-Newton, multiplicative updates, and all-at-once methods: it is the dominant factor-update computation itself (Phan et al., 2012).

The same centrality persists in generalized CP (GCP), where the objective is separable over tensor entries and may use Bernoulli, Huber, or other losses. If

$F_0, F_1, \ldots, F_{N-1}$9

then the exact gradient with respect to $$\mathcal{X} \approx [\![\lambda; F_0, F_1, \ldots, F_{N-1}]\!] = \sum_{r=0}^{R-1} \lambda_r \, f^{(0)}_r \circ f^{(1)}_r \circ \cdots \circ f^{(N-1)}_r.$$0 is

$$\mathcal{X} \approx [\![\lambda; F_0, F_1, \ldots, F_{N-1}]\!] = \sum_{r=0}^{R-1} \lambda_r \, f^{(0)}_r \circ f^{(1)}_r \circ \cdots \circ f^{(N-1)}_r.$$1

with $$\mathcal{X} \approx [\![\lambda; F_0, F_1, \ldots, F_{N-1}]\!] = \sum_{r=0}^{R-1} \lambda_r \, f^{(0)}_r \circ f^{(1)}_r \circ \cdots \circ f^{(N-1)}_r.$$2 the Khatri-Rao product of all factor matrices except $$\mathcal{X} \approx [\![\lambda; F_0, F_1, \ldots, F_{N-1}]\!] = \sum_{r=0}^{R-1} \lambda_r \, f^{(0)}_r \circ f^{(1)}_r \circ \cdots \circ f^{(N-1)}_r.$$3. The key difficulty is that the derivative tensor $$\mathcal{X} \approx [\![\lambda; F_0, F_1, \ldots, F_{N-1}]\!] = \sum_{r=0}^{R-1} \lambda_r \, f^{(0)}_r \circ f^{(1)}_r \circ \cdots \circ f^{(N-1)}_r.$$4 is dense almost everywhere even if $$\mathcal{X} \approx [\![\lambda; F_0, F_1, \ldots, F_{N-1}]\!] = \sum_{r=0}^{R-1} \lambda_r \, f^{(0)}_r \circ f^{(1)}_r \circ \cdots \circ f^{(N-1)}_r.$$5 is sparse. The stochastic construction replaces $$\mathcal{X} \approx [\![\lambda; F_0, F_1, \ldots, F_{N-1}]\!] = \sum_{r=0}^{R-1} \lambda_r \, f^{(0)}_r \circ f^{(1)}_r \circ \cdots \circ f^{(N-1)}_r.$$6 by a sparse unbiased estimator $$\mathcal{X} \approx [\![\lambda; F_0, F_1, \ldots, F_{N-1}]\!] = \sum_{r=0}^{R-1} \lambda_r \, f^{(0)}_r \circ f^{(1)}_r \circ \cdots \circ f^{(N-1)}_r.$$7 satisfying

$$\mathcal{X} \approx [\![\lambda; F_0, F_1, \ldots, F_{N-1}]\!] = \sum_{r=0}^{R-1} \lambda_r \, f^{(0)}_r \circ f^{(1)}_r \circ \cdots \circ f^{(N-1)}_r.$$8

so that

$$\mathcal{X} \approx [\![\lambda; F_0, F_1, \ldots, F_{N-1}]\!] = \sum_{r=0}^{R-1} \lambda_r \, f^{(0)}_r \circ f^{(1)}_r \circ \cdots \circ f^{(N-1)}_r.$$9

This converts a dense gradient evaluation of cost $n$0 into a sparse MTTKRP of cost $n$1 with storage $n$2 (Kolda et al., 2019).

In distributed nonnegative CP decomposition, MTTKRP similarly supplies the right-hand side for the local nonnegative least-squares subproblems. For mode $n$3,

$n$4

After $n$5 is formed, the row-wise NNLS updates are local, but the expensive global step remains MTTKRP. This is precisely why dimension trees are introduced in dense NNCP: to reuse partial contractions across multiple MTTKRPs within one alternating cycle (Ballard et al., 2018).

3. Bottleneck structure, irregularity, and communication

The sparse kernel is expensive for reasons that differ sharply from dense matrix multiplication. In sparse CP-ALS, each nonzero touches several factor-matrix rows, updates are scattered across output rows, the same output row may be updated by many nonzeros, and synchronization is needed to prevent races. The result is a kernel described as memory-bound, irregular, and synchronization-heavy rather than a dense linear algebra kernel (Wijeratne et al., 20 Jul 2025).

This characterization is echoed in FPGA-oriented analyses, which emphasize that each nonzero requires multiple memory reads for tensor coordinates, tensor values, factor-matrix rows, and read-modify-write access to the output row or fiber. The access stream contains many secondary cache misses, multiple outstanding misses can occur to the same cache line, and the dominant cost is external memory latency and bandwidth rather than multiply-add throughput (Wijeratne et al., 2021).

At very large scale, additional system effects become first-order. Billion-scale sparse tensors can exceed the global memory of a single GPU, especially when implementations require multiple tensor copies for different mode orderings. Prior single-GPU approaches then encounter memory-capacity limits, irregular access imbalance, communication overhead for sharing factor matrices across modes, and synchronization overhead when multiple workers update the same output rows (Wijeratne et al., 20 Jul 2025).

For dense tensors, the dominant issue is formalized as data movement. Communication lower bounds for MTTKRP show that the structure of the computation allows less communication than the straightforward approach of casting the computation as a matrix multiplication operation. Sequential and parallel algorithms can attain these lower bounds, establishing that tensor-aware MTTKRP is not simply a reformulation of GEMM, but a distinct communication problem with exploitable structure (Ballard et al., 2017).

A common misconception is therefore that MTTKRP is reducible, without penalty, to explicit Khatri-Rao formation followed by ordinary matrix multiplication. The literature repeatedly rejects that view. Dense matrix-free implementations in GenTen contrast with matrix-based methods used by Tensor Toolbox and TensorLy, while communication-theoretic analyses show that tensor structure can reduce data movement below that of the corresponding matrix product (Kosmacher et al., 16 Oct 2025, Ballard et al., 2017).

4. Algorithmic reformulations and computational reuse

Several major lines of work seek to accelerate MTTKRP by reorganizing the computation rather than changing the objective. One dense strategy avoids explicit interior-mode unfoldings and instead rewrites the kernel as recursive tensor-vector contractions. The method chooses a left-to-right or right-to-left contraction order based on tensor dimensions and reuses intermediate contractions across adjacent modes. Reported speedups are about $n$6 for 3-D tensors, about $n$7 for 4-D tensors, often $n$8 for higher-order tensors, and $n$9 in a 7-D example (Phan et al., 2012).

A related shared-memory line casts nearly all dense MTTKRP arithmetic as matrix operations while avoiding tensor-entry reordering in memory. The 1-step and 2-step algorithms exploit BLAS-friendly block structure and partial MTTKRP plus multi-TTV decompositions. In a neuroimaging CP decomposition, the implementation achieved up to $$X_{(n)} \left( F_{N-1} \odot \cdots \odot F_{n+1} \odot F_{n-1} \odot \cdots \odot F_0 \right),$$0 speedup over existing parallel software (Hayashi et al., 2017).

Dimension trees attack repeated work across modes in dense alternating optimization. In distributed-memory nonnegative CP decomposition, the tree organizes the tensor at the root, partial MTTKRPs at internal branches, and final mode-specific outputs at the leaves. The paper states that this can reduce computation by roughly a factor of $$X_{(n)} \left( F_{N-1} \odot \cdots \odot F_{n+1} \odot F_{n-1} \odot \cdots \odot F_0 \right),$$1 for dense $$X_{(n)} \left( F_{N-1} \odot \cdots \odot F_{n+1} \odot F_{n-1} \odot \cdots \odot F_0 \right),$$2-way tensors in the best case, and reports the dimension-tree version as $$X_{(n)} \left( F_{N-1} \odot \cdots \odot F_{n+1} \odot F_{n-1} \odot \cdots \odot F_0 \right),$$3 faster than the version without dimension trees on a 5D synthetic tensor (Ballard et al., 2018).

Compiler-oriented workspaces provide a different kind of reformulation, especially for sparse kernels. For the 3D expression

$$X_{(n)} \left( F_{N-1} \odot \cdots \odot F_{n+1} \odot F_{n-1} \odot \cdots \odot F_0 \right),$$4

workspace transformations hoist loop-invariant subexpressions such as $$X_{(n)} \left( F_{N-1} \odot \cdots \odot F_{n+1} \odot F_{n-1} \odot \cdots \odot F_0 \right),$$5 into temporaries and avoid expensive insertion into sparse results. This led to a $$X_{(n)} \left( F_{N-1} \odot \cdots \odot F_{n+1} \odot F_{n-1} \odot \cdots \odot F_0 \right),$$6 improvement for MTTKRP on the nell-2 data set, and also enabled sparse-output MTTKRP; at density $$X_{(n)} \left( F_{N-1} \odot \cdots \odot F_{n+1} \odot F_{n-1} \odot \cdots \odot F_0 \right),$$7, the sparse-output variant reported $$X_{(n)} \left( F_{N-1} \odot \cdots \odot F_{n+1} \odot F_{n-1} \odot \cdots \odot F_0 \right),$$8 to $$X_{(n)} \left( F_{N-1} \odot \cdots \odot F_{n+1} \odot F_{n-1} \odot \cdots \odot F_0 \right),$$9 speedups depending on the tensor (Kjolstad et al., 2018).

The recent dense matrix-free line in GenTen pushes the same idea further by eliminating explicit Khatri-Rao materialization altogether. The matrix-free memory cost grows with rank like the sum of tensor mode sizes, $X_{(n)}$0, whereas matrix-based methods grow like the Khatri-Rao expansion. For the largest reported problem, a rank-2000 MTTKRP, the matrix-free memory cost was about $X_{(n)}$1 of the matrix-based methods, a $X_{(n)}$2 reduction. The optimized matrix-free implementation reported a $X_{(n)}$3 single-device speedup on an Intel 8480+ CPU and an $X_{(n)}$4 speedup on an NVIDIA H100 GPU over baseline matrix-free implementations, and enabled a rank-2000 decomposition on a single H100 instead of six H100s for the matrix-based variant (Kosmacher et al., 16 Oct 2025).

5. Parallel and architectural realizations

The implementation literature is unusually heterogeneous because the kernel’s limiting factor changes with memory hierarchy, synchronization model, and device scale. The following systems illustrate distinct design points.

Context	Representative mechanism	Reported outcome
Multi-core CPU	Dynasor uses FLYCOO, dynamic tensor remapping, lock-free concurrent processing of independent partitions, and greedy super-shard scheduling	$X_{(n)}$5 to $X_{(n)}$6 speedup in total execution time across all modes versus state-of-the-art CPU implementations (Wijeratne et al., 2023)
Single GPU	Dynamic tensor remapping groups dependent nonzeros into independent thread-block partitions and avoids global atomics across thread blocks	$X_{(n)}$7, $X_{(n)}$8, and $X_{(n)}$9 geometric mean speedup versus BLCO, MM-CSF, and ParTI-GPU (Wijeratne et al., 2024)
Multi-GPU	AMPED partitions nonzeros by output-mode index, load-balances inter-device partitions greedily, and uses ring-based GPUDirect P2P communication	$n$0 geometric mean speedup in total execution time over state-of-the-art GPU baselines using 4 GPUs on billion-scale tensors (Wijeratne et al., 20 Jul 2025)
FPGA memory system	Reconfigurable Local Memory Blocks combine non-blocking cache, DMA engine, and Request Reductor	$n$1 speedup over commercial memory controller IPs, $n$2 over cache-only, and $n$3 over DMA-only (Wijeratne et al., 2021)
PIM and photonic memory	PRISM maps spMTTKRP to UPMEM PIM with fixed-point arithmetic and heterogeneous CPU collaboration; photonic pSRAM maps Hadamard-scale-accumulate primitives onto WDM-enabled in-memory compute	Up to $n$4 PIM-only and $n$5 heterogeneous speedup over state-of-the-art CPU implementations; a separate predictive pSRAM model reports 17 PetaOps sustained MTTKRP performance (Pacheco et al., 28 May 2026, Wijeratne et al., 23 Mar 2025)

These systems share several recurrent principles. One is output-row ownership: if all nonzeros contributing to the same output row are kept in one partition, inter-worker races can be removed or confined. AMPED uses output-mode-index partitioning to ensure that no two GPUs update the same output row during a mode, while CPU and GPU dynamic-remapping schemes pursue the same effect at thread or thread-block granularity (Wijeratne et al., 20 Jul 2025, Wijeratne et al., 2023).

Another is explicit treatment of MTTKRP as a memory-system problem. FPGA work does not primarily modify arithmetic; it redesigns the memory substrate through cache, DMA, request coalescing, or even optical SRAM. This suggests that, for sparse MTTKRP, external-memory organization is often more consequential than raw compute throughput (Wijeratne et al., 2021, Wijeratne et al., 2022).

The system literature also qualifies simple scaling narratives. AMPED reports nearly linear compute scaling with the number of GPUs on large tensors, yet also notes that the single-GPU FLYCOO-GPU can outperform it on the smallest tensor, Twitch, because communication overhead dominates there. This makes scale sensitivity an intrinsic property of the kernel rather than an implementation accident (Wijeratne et al., 20 Jul 2025).

6. Sampling, portability, and current directions

A major recent direction is to accelerate MTTKRP-adjacent least-squares solves by sampling the Khatri-Rao structure instead of redesigning the core accumulation kernel. Exact leverage-score sampling from Khatri-Rao products makes it possible to sketch the CP-ALS least-squares subproblems without materializing the exponentially tall design matrix. The resulting STS-CP method is reported to achieve $n$6 of the fit of non-randomized ALS faster than CP-ARLS-LEV, with about $n$7 to $n$8 speedup on Amazon and Reddit; on Enron, it used about 65,000 samples per solve beyond rank 50 while CP-ARLS-LEV could not achieve the same accuracy even with many more samples (Bharadwaj et al., 2023).

Sampling also appears one level earlier, at the derivative-tensor level. In stochastic GCP, the algorithm does not approximate the algebraic form of the gradient; it approximates the dense derivative tensor $n$9 by a sparse unbiased $\odot$0, after which the exact same MTTKRP structure applies. This is a notable conceptual point: sampling and MTTKRP are complementary rather than competing strategies (Kolda et al., 2019).

Performance portability remains unsettled. In the PASTA MTTKRP benchmark, a Kokkos port achieved geometric mean speedups of about $\odot$1 to $\odot$2 over the PASTA baseline on CPU systems, but showed a $\odot$3 slowdown on NVIDIA A100 relative to the PASTA CUDA baseline. The same study attributed low fractions of theoretical peak memory bandwidth not to a failure of Kokkos, but to the fact that MTTKRP is bottlenecked by memory load/store behavior and is latency-bound. Separately, the analysis of the related $\odot$4 kernel in SparTen concluded that atomic operations were not a critical bottleneck while higher cache reuse could provide a non-trivial performance improvement (Anderson et al., 2023).

This broader record suggests that MTTKRP should be understood less as a single algorithm than as a family of contraction/update kernels whose optimal formulation depends on tensor density, tensor order, rank, memory hierarchy, synchronization scope, and the cost of explicit Khatri-Rao formation. Dense and sparse cases diverge sharply; single-device and distributed settings diverge again; and emerging hardware proposals increasingly treat MTTKRP as a memory-centric workload rather than a purely arithmetic one. Across these settings, however, one feature remains constant: MTTKRP is the operation that couples tensor data to the current factor matrices, and therefore the place where tensor decomposition performance is won or lost.