Extension of NMF to 3D AFAIK, Non-Negative Matrix Factorization (NMF) is the procedure of looking for matrices $A$ and $B$ such that

$$Data_{ik} = \sum_j A_{ij} B_{jk}$$

My data matrix is in fact 3D. I would like to fit the following model to my data

$$Data_{ikl} = \sum_j A_{ij} B_{jk} C_{jl}$$

It would be great to know if this model has a name and is already implemented somewhere (preferably Python)
Are you looking for something like this?

"Introduction to Tensor Decompositions and their
Applications in Machine Learning"
https://arxiv.org/pdf/1711.10781.pdf

I'm not sure if there are better papers for it, but from what I can tell tensor decomposition seems like what you're describing.