AbstractModern instrumentation in chemistry routinely generates two‐dimensional (second‐order) arrays of data. Considering that most analyses need to compare several samples, the analyst ends up with a three‐dimensional (third‐order) array which is difficult to visualize or interpret with the conventional statistical tools.Some of these data arrays follow the so‐calledtrilinearmodel,\documentclass{article}\pagestyle{empty}\begin{document}$$ {\rm R}_{ijk} = \sum\limits_{r = 1}^N {{\rm X}_{ir} {\rm Y}_{jr} {\rm Z}_{kr} + {\rm Error}_{ijk} } $$\end{document}These trilinear arrays of data are known to have unique factor analysis decompositions which correspond to the true physical factors that form the data, i.e. given the array ℝ, a unique solution can be found in many cases for each orderX,YandZ. This is in contrast to the well‐known second‐order bilinear data factor analysis, where the abstract solutions obtained are not unique and at best cannot be easily compared with the underlying physical factors owing to a rotational ambiguity.Trilinear decompositions have had the disadvantage, however, that a non‐linear optimization with many parameters is necessary to reach a least‐squares solution. This paper will introduce a method for reducing the problem to a rectangular generalized eigenvalue–eigenvector equation where the eigenvectors are the contravariant form (pseudo‐inverse) of the actual factors. It is shown that the method works well when the factors are linearly independent in at least two orders (e.g.XirandYjrare full rank matrices).Finally, it is shown how trilinear decompositions relate to multicomponent calibration, curve resolutio