In this article, we address the problem of tensor factorization subject to certain constraints. We focus on the canonical polyadic decomposition, also known as parallel factor analysis. The interest of this multilinear decomposition coupled with 3D fluorescence spectroscopy is now well established in the fields of environmental data analysis, biochemistry, and chemistry. When real experimental data (possibly corrupted by noise) are processed, the actual rank of the “observed” tensor is generally unknown. Moreover, when the amount of data is very large, this inverse problem may become numerically ill‐posed and consequently hard to solve. The use of proper constraints reflecting some a priori knowledge about the latent (or hidden) tracked variables and/or additional information through the addition of penalty functions can prove very helpful in estimating more relevant components rather than totally arbitrary ones. The counterpart is that the cost functions that have to be considered can be nonconvex and sometimes even nondifferentiable, making their optimization more difficult and leading to a higher computing time and a slower convergence speed. Block alternating proximal approaches offer a rigorous and flexible framework to properly address that problem since they are applicable to a large class of cost functions while remaining quite easy to implement. Here, we suggest a new block coordinate variable metric forward‐backward method, which can be seen as a special case of majorize‐minimize approaches to derive a new penalized nonnegative third‐order canonical polyadic decomposition algorithm. Its interest, efficiency, robustness, and flexibility are illustrated thanks to computer simulations performed on both simulated and real experimental 3D fluorescence spectroscopy data.