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)