Parallel Factor Analysis

PARAFAC is a decomposition of a multiway array into a sum of rank-one tensors. Harshman introduced the procedure in 1970 as Parallel Factor Analysis, an explanatory multimodal factor model built from Cattell's principle of rotation to proportional profiles. Independently in the same year, Carroll and Chang introduced CANDECOMP, a canonical decomposition for N-way generalizations of the Eckart-Young model in individual-differences scaling. Kolda and Bader, following Kiers, refer to the common decomposition as CANDECOMP/PARAFAC, abbreviated CP. This page uses PARAFAC for the chemometric name and CP for the same mathematical model.

The V1 treatment is classical three-way CP fitted by unconstrained alternating least squares. Higher-order arrays, missing values, nonnegativity, PARAFAC2, and line search are outside the educational algorithm. TensorLy and the N-way Toolbox are software examples, not mathematical proofs.

Kolda and Bader define an N-way, or Nth-order, tensor as a multidimensional array. A first-order tensor is a vector, a second-order tensor is a matrix, and order three or higher is a higher-order tensor. The dimensions are also called modes or ways. A three-way array is written

with entries $x_{ijk}$. Not every collection of matrices is an appropriate PARAFAC data structure. Bro emphasizes that the array must support a multilinear model. Unfolding a cube into a matrix and running PCA is a different model, not PARAFAC.

One chemometric three-way structure is fluorescence excitation-emission data arranged as samples by excitation wavelengths by emission wavelengths. Stedmon and Bro describe excitation-emission matrices (EEMs) as well suited to multiway analysis, including PARAFAC. That example is a data layout, not a claim that every PARAFAC component is a pure fluorophore.

For $\mathcal{X}\in {\mathbb{R}}^{I\times J\times K}$, Kolda and Bader write the CP approximation as a sum of $R$ rank-one tensors

Bro writes the same trilinear form with an explicit residual for measured data,

The residual tensor $\mathcal{E}$ with entries $e_{ijk}$ is the difference between the data and the reconstruction. It is not automatically noise. It can contain noise, unmodelled components, and departures from multilinearity.

Kolda and Bader also use a normalized convention

in which the columns of the factor matrices have unit Euclidean length and the scale of each rank-one term is stored in ${\lambda }_r$. The educational algorithm follows that convention. It represents the same family of reconstructions as unnormalized factor matrices whose column scales multiply to the same product.

Kolda and Bader define a third-order tensor as rank one when it is an outer product of three vectors,

with entries $x_{ijk}=a_i{b}_j{c}_k$. The symbol $\circ$ is the vector outer product. It is not a dot product, not a Kronecker product $\otimes$, and not a Khatri-Rao product $\odot$.

Tensor rank is the smallest number of rank-one tensors that sum exactly to $\mathcal{X}$. Kolda and Bader stress that this rank does not behave like matrix rank: there is no Eckart-Young truncation of a CP decomposition, determining tensor rank is NP-hard in general, and a best rank-$R$ approximation need not exist. Typical rank and border rank are advanced remarks only. For noisy chemometric data, the fitted $R$ is a model-selection choice, not a computed matrix rank.

The factor matrices collect the rank-one vectors as columns:

Columns with the same index $r$ belong to one component. Bro notes that three-way practice often does not distinguish scores from loadings because the three modes are treated equally in the numerical model. This page does not call $\mathbf{A}$ scores universally, and it does not call $\mathbf{B}$ or $\mathbf{C}$ spectra universally. Interpretation follows what each mode measures.

In fluorescence EEM data, one factor matrix can describe sample-dependent component contributions, one can describe excitation-mode profiles, and one can describe emission-mode profiles. Stedmon and Bro present PARAFAC on DOM fluorescence as a tutorial for that layout. Sample factors are not automatically absolute concentrations. Scale is a modelling convention unless the analytical experiment identifies it.

PARAFAC factor columns are not required to be orthogonal. Kolda and Bader contrast this with SVD, whose uniqueness in the matrix case depends on orthogonality. Orthogonality is an optional constraint in some software. It is not part of the unconstrained CP model taught here.

Mode-n matricization arranges the mode-n fibres as the columns of a matrix ${\mathbf{X}}_{(n)}$. This page uses the Kolda and Bader convention, including their worked 3 by 4 by 2 example. Different papers and packages permute those columns differently. Kolda and Bader state that the permutation is not important provided it is consistent across related calculations. TensorLy documents a different unfolding chosen for C-order arrays. The educational Python and MATLAB code do not use the TensorLy convention.

For $\mathcal{X}\in {\mathbb{R}}^{I\times J\times K}$:

• ${\mathbf{X}}_{(1)}\in {\mathbb{R}}^{I\times JK}$ has columns ${\mathbf{x}}_{:jk}$ with $j$ cycling faster than $k$.
• ${\mathbf{X}}_{(2)}\in {\mathbb{R}}^{J\times IK}$ has columns ${\mathbf{x}}_{i:k}$ with $i$ cycling faster than $k$.
• ${\mathbf{X}}_{(3)}\in {\mathbb{R}}^{K\times IJ}$ has columns ${\mathbf{x}}_{ij:}$ with $i$ cycling faster than $j$.

The Khatri-Rao product $\mathbf{A}\odot \mathbf{B}$ of $\mathbf{A}\in {\mathbb{R}}^{I\times R}$ and $\mathbf{B}\in {\mathbb{R}}^{J\times R}$ is the matching columnwise Kronecker product

with ${\mathbf{b}}_r$ inner in ${\mathbf{a}}_r\otimes {\mathbf{b}}_r$. It is not the full Kronecker product $\otimes$ of the two matrices, and it is not a Hadamard (entrywise) product $\ast$.

Under this unfolding and this Khatri-Rao order, Kolda and Bader give

Dimension check for mode 1: $\mathbf{C}\odot \mathbf{B}$ is $(KJ)\times R$, so $\mathbf{A}(\mathbf{C}\odot \mathbf{B}{)}^{\top }$ is $I\times JK$, matching ${\mathbf{X}}_{(1)}$. Modes 2 and 3 follow analogously. Reconstruction in the educational code uses an independent outer-product sum, not the unfolding helpers, so a consistent identity can be tested rather than assumed.

For fixed $R$, Kolda and Bader state the least-squares CP problem as minimizing $\Vert \mathcal{X}-\hat{\mathcal{X}}{\Vert }_F$ over the factor representation of $\hat{\mathcal{X}}$. The joint problem is not a single linear least-squares problem and is not convex. Alternating least squares, proposed in the original Carroll-Chang and Harshman papers and presented by Kolda and Bader as the workhorse CP algorithm, fixes all but one factor matrix and solves the remaining linear least-squares subproblem.

With $\mathbf{B}$ and $\mathbf{C}$ fixed,

The educational code solves that subproblem with a linear least-squares solver, equivalent to

Kolda and Bader also give the Gram form

where $\ast$ is the Hadamard product. They note that reducing the pseudoinverse to an $R\times R$ matrix can be numerically ill conditioned. The educational implementation therefore does not use an explicit inverse, and it does not use that reduced Gram formula.

The mode-2 and mode-3 updates are the analogous least-squares problems on ${\mathbf{X}}_{(2)}$ and ${\mathbf{X}}_{(3)}$, using $\mathbf{C}\odot \mathbf{A}$ and $\mathbf{B}\odot \mathbf{A}$. After each mode update, columns are normalized and the column norms replace $\mathbf{\lambda }$, as in Kolda and Bader Figure 3.3. All $R$ components are estimated together. A rank-one fit followed by deflation is not standard PARAFAC. Kolda and Bader state that a best rank-k tensor approximation does not generally arise by taking leading CP terms one at a time.

ALS is iterative. Kolda and Bader state that it is not guaranteed to converge to a global minimum or even a stationary point, only to a solution where the objective ceases to decrease, and that the result can depend heavily on the starting guess. Bro likewise notes that ALS improves (or does not worsen) the fit at each iteration, and that convergence to the global least-squares solution is what one hopes for in well-behaved problems, not a theorem. TensorLy documentation that speaks of a reconstruction tolerance as indicating a global minimum is software wording. It is not used here as a mathematical guarantee.

Random starts and SVD-based starts from mode unfoldings are both described by Kolda and Bader. Harshman and Lundy, as reported by Bro, advocate several random starts. Agreement of several starts reduces, but does not remove, the chance of an unfortunate local solution. The lowest reconstruction error among tested starts is the best solution among those starts. It is not proved to be the global minimizer.

Classical three-way PARAFAC on this page is the Kolda-Bader CP model, the three matricized identities, unconstrained ALS subproblems, column normalization into $\mathbf{\lambda }$, and Kruskal's sufficient uniqueness condition.

Equation (2) is Bro's residual form with unnormalized factors; it is the same reconstruction as (1) after absorbing ${\lambda }_r$ into one of the three vectors. Equation (8) is a sufficient uniqueness condition, not a noisy-data diagnostic. Equation (9) is the relative reconstruction error used in the educational code. It is not labelled $R^2$ and it is not a package-specific fit statistic.

The public algorithm matches Kolda and Bader Figure 3.3 restricted to three modes, with an explicit relative-error change stop. That stop belongs to their class of little-or-no-objective-improvement rules. It is an educational choice, not a universal PARAFAC constant, and not the TensorLy default presented as science.

Scale is not separately identified. Kolda and Bader state that

represents the same rank-one term whenever ${\alpha }_r{\beta }_r{\gamma }_r=1$. Column normalization therefore stores a convention, not extra chemical information. Unconstrained real factors can also redistribute sign among the three vectors of a component while keeping that product equal to one.

Component order is arbitrary. Any simultaneous permutation of corresponding columns of $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{C}$ leaves the reconstruction unchanged. Two solutions must be matched before column 1 is compared with column 1.

Essential uniqueness means uniqueness up to those scale and permutation indeterminacies. It does not mean that PARAFAC always recovers true components, and it does not mean that there is no ambiguity. Bro contrasts this with the ordinary rotational freedom of unconstrained bilinear factor models: if a PARAFAC decomposition is essentially unique, it does not possess that arbitrary rotational indeterminacy. That statement is conditional on identifiability conditions and on an appropriate trilinear model. It is not a claim that PARAFAC has no ambiguity.

Kruskal k-rank $k_{\mathbf{A}}$ is the largest integer $k$ such that every subset of $k$ columns of $\mathbf{A}$ is linearly independent. It is not ordinary matrix rank. Kolda and Bader state Kruskal's sufficient condition for essential uniqueness of a three-way rank-$R$ CP decomposition as

This page uses that Kolda-Bader statement. Bro 1997 writes a related k-rank inequality with a strict greater-than sign. The two forms are not silently merged. The condition is sufficient, not necessary in general. Ten Berge and Sidiropoulos, as reported by Kolda and Bader, showed necessity for $R=2$ and $R=3$ but not for $R>3$. The inequality is not a diagnostic for noisy fitted data. Sidiropoulos and Bro extend uniqueness to N-way arrays; that extension is not derived here.

Fitted $R$ is not matrix rank of an unfolding, not a 95% variance cutoff, and not a single singular-value gap. Bro states that deciding the rank of a PARAFAC model is difficult, and that successive deflation as in PCA is not the PARAFAC construction. Over-factoring can split true factors across correlated components rather than merely adding a noise component.

Evidence can include reconstruction, residual structure, stability across starts, core consistency, split-half checks in some fluorescence studies, and external chemical knowledge. None of these is a universal oracle.

Bro and Kiers introduced the core consistency diagnostic, CORCONDIA, for assessing whether a fitted PARAFAC model remains compatible with the restricted superdiagonal interaction structure expected from PARAFAC when compared with a more flexible Tucker3 representation. A PARAFAC model is called appropriate, in their wording, if adding other combinations of the same components does not improve the fit considerably. They propose choosing the largest model that is still sufficiently appropriate. They also state that the theoretical understanding of CORCONDIA is not yet complete. This page does not implement CORCONDIA, does not reproduce an equation from a secondary source, and does not teach a universal threshold such as 90% correct or 50% invalid.

Split-half analysis appears in Bro, following Harshman and Lundy, and in Stedmon and Bro for fluorescence EEMs. It is an application strategy in that literature: comparable loadings in the unsplit modes, after permutation and scale matching, support a candidate $R$. It is not the universal validation procedure for all PARAFAC applications.

Kolda and Bader review degeneracy: a tensor may be approximated arbitrarily well by a factorization of lower rank, in which case a best rank-k approximation need not exist. Practical symptoms discussed by Bro and by Kolda and Bader include very large component magnitudes of opposite sign, nearly proportional factors, and unstable or nonconvergent fits. Bro monitors a triple cosine between component triples as a degeneracy indication. Not every slow ALS run is degeneracy. Causes discussed by Bro include extracting too many components, poor preprocessing, and an inappropriate (non-trilinear) model.

Kolda and Bader list stopping on little objective improvement, little factor change, an objective near zero, or a maximum iteration count. TensorLy currently exposes n_iter_max, tol, and cvg_criterion with options abs_rec_error and rec_error. Those names are TensorLy behaviour. TensorLy documentation that treats a reconstruction-error tolerance as indicating a global minimum is software wording, not a mathematical guarantee. The educational stop is the absolute change in relative reconstruction error. Software defaults are not scientific recommendations.

A small reconstruction error does not by itself demonstrate correct $R$, essential uniqueness, chemical validity, correct component identities, absence of degeneracy, or robustness. Fit is not model validity.

Bro presents PARAFAC as a multiway method originating in psychometrics and used in chemometrics for crossed multiway measurements, including fluorescence emission spectra at several excitation wavelengths for several samples. Under an appropriate multilinear measurement model and sufficient identifiability conditions, PARAFAC can yield factors with direct chemical interpretability. Bro discusses recovery of pure spectra in suitable multiway spectroscopic conditions. That is not a universal statement. PARAFAC does not assign chemical names. A resolved loading is not a molecular identification. Stedmon and Bro warn that interpretation of DOM fluorescence components can be more complex than assigning each factor to a pure fluorophore.

Bro emphasizes that scaling and centering multiway data are not as straightforward as in two-way analysis. Centering is across a mode; scaling is of slabs, not ordinary two-way autoscaling of mixed unfoldings. This page does not teach a universal "always autoscale" or "always mean-center each mode" rule.

PCA operates on a two-way matrix and produces orthogonal components ordered by captured variance under the PCA/SVD construction. Bro calls PARAFAC a generalization of PCA to higher order arrays and immediately stresses differences: no ordinary rotation problem when the model is essentially unique, components are not estimated successively, and preprocessing is not the two-way default. PARAFAC is not PCA for tensors.

SVD decomposes a matrix. The Eckart-Young optimal truncation of leading singular components does not carry over to CP. Kolda and Bader state that explicitly.

Tucker uses factor matrices plus a generally dense core. Kiers, as used by Bro, and Kolda and Bader, treat CP/PARAFAC as a restricted Tucker structure with a superdiagonal core and matching component counts across modes. Tucker is not PARAFAC. Not every tensor decomposition is PARAFAC.

MCR-ALS typically models a two-way bilinear matrix $\mathbf{D}=\mathbf{C}{\mathbf{S}}^{\top }+\mathbf{E}$ and uses constraints to obtain interpretable profiles. PARAFAC uses a multiway multilinear structure: one vector in each mode linked by an outer product. Neither method is universally superior. Rotational ambiguity is the two-way bilinear issue; essential uniqueness, when it holds, is the three-way CP property discussed above.

The first Python listing is unconstrained three-way CP ALS under the Kolda and Bader unfolding and Khatri-Rao order. Updates use numpy.linalg.lstsq. Reconstruction uses numpy.einsum independently of unfolding. After each mode update, column norms replace $\mathbf{\lambda }$. The example tolerance tol=1e-8 is an educational relative-error change threshold, not a universal default. Missing values, nonnegativity, CORCONDIA, and line search are not implemented.

The second listing calls current TensorLy parafac and reconstructs with cp_to_tensor. TensorLy returns a CPTensor of weights and factor matrices. TensorLy unfolding differs from Kolda and Bader. TensorLy init, tol, and cvg_criterion are package options. They do not prove a global minimum.

MATLAB implements the same unfolding, the same Khatri-Rao helper with kron column by column, and transpose-backslash least squares. Base MATLAB has no khatrirao. The N-way Toolbox of Andersson and Bro (current official University of Copenhagen Chemometrics Research release 3.60, September 2026) provides PARAFAC, Tucker, N-PLS, constraints, multiway centering and scaling, cross-validation, and core consistency. Its parafac calling convention is not copied here from memory. Tensor Toolbox cp_als is a separate independent MATLAB CP ALS implementation. It is not used in the educational listing.

import numpy as np  def _validate_tensor(X, label="X"):    X = np.asarray(X, dtype=float)    if X.ndim != 3:        raise ValueError("%s must be a 3-way tensor." % label)    if min(X.shape) < 1:        raise ValueError("%s must have positive dimensions." % label)    if not np.all(np.isfinite(X)):        raise ValueError("%s must contain only finite values." % label)    return X  def _validate_factor(F, n_rows, rank, label):    F = np.asarray(F, dtype=float)    if F.ndim != 2:        raise ValueError("%s must be a 2D array." % label)    if F.shape != (n_rows, rank):        raise ValueError("%s must have shape (%d, %d)." % (label, n_rows, rank))    if not np.all(np.isfinite(F)):        raise ValueError("%s must contain only finite values." % label)    return F  def unfold(X, mode):    """Kolda and Bader (2009) mode-n unfolding. mode is 1, 2, or 3."""    X = _validate_tensor(X)    I, J, K = X.shape    if mode == 1:        return np.column_stack(            [X[:, j, k] for k in range(K) for j in range(J)]        )    if mode == 2:        return np.column_stack(            [X[i, :, k] for k in range(K) for i in range(I)]        )    if mode == 3:        return np.column_stack(            [X[i, j, :] for j in range(J) for i in range(I)]        )    raise ValueError("mode must be 1, 2, or 3.")  def khatri_rao(A, B):    """Matching columnwise Kronecker product A odot B (Kolda and Bader)."""    A = np.asarray(A, dtype=float)    B = np.asarray(B, dtype=float)    if A.ndim != 2 or B.ndim != 2:        raise ValueError("Khatri-Rao factors must be 2D arrays.")    if A.shape[1] != B.shape[1]:        raise ValueError("Khatri-Rao factors must share the column count R.")    rank = A.shape[1]    out = np.empty((A.shape[0] * B.shape[0], rank), dtype=float)    for r in range(rank):        out[:, r] = np.kron(A[:, r], B[:, r])    return out  def reconstruct(weights, A, B, C):    weights = np.asarray(weights, dtype=float).reshape(-1)    A = np.asarray(A, dtype=float)    B = np.asarray(B, dtype=float)    C = np.asarray(C, dtype=float)    if not (A.shape[1] == B.shape[1] == C.shape[1] == weights.shape[0]):        raise ValueError("weights and factor matrices must share R.")    return np.einsum("r,ir,jr,kr->ijk", weights, A, B, C)  def relative_reconstruction_error(X, Xhat):    X = _validate_tensor(X)    Xhat = np.asarray(Xhat, dtype=float)    if Xhat.shape != X.shape:        raise ValueError("Xhat must have the same shape as X.")    den = np.linalg.norm(X)    if den == 0:        raise ValueError("X must have a nonzero Frobenius norm.")    return float(np.linalg.norm(X - Xhat) / den)  def normalize_columns(F):    F = np.asarray(F, dtype=float)    norms = np.linalg.norm(F, axis=0)    norms = np.where(norms == 0.0, 1.0, norms)    return F / norms, norms  def update_A(X, B, C):    X = _validate_tensor(X)    rank = B.shape[1]    B = _validate_factor(B, X.shape[1], rank, "B")    C = _validate_factor(C, X.shape[2], rank, "C")    Z = khatri_rao(C, B)    AT, *_ = np.linalg.lstsq(Z, unfold(X, 1).T, rcond=None)    return AT.T  def update_B(X, A, C):    X = _validate_tensor(X)    rank = A.shape[1]    A = _validate_factor(A, X.shape[0], rank, "A")    C = _validate_factor(C, X.shape[2], rank, "C")    Z = khatri_rao(C, A)    BT, *_ = np.linalg.lstsq(Z, unfold(X, 2).T, rcond=None)    return BT.T  def update_C(X, A, B):    X = _validate_tensor(X)    rank = A.shape[1]    A = _validate_factor(A, X.shape[0], rank, "A")    B = _validate_factor(B, X.shape[1], rank, "B")    Z = khatri_rao(B, A)    CT, *_ = np.linalg.lstsq(Z, unfold(X, 3).T, rcond=None)    return CT.T  def parafac_als(    X,    rank,    max_iter=100,    tol=1e-8,    random_state=None,    A0=None,    B0=None,    C0=None,):    X = _validate_tensor(X)    rank = int(rank)    if rank < 1:        raise ValueError("rank R must be at least 1.")    I, J, K = X.shape    if A0 is None or B0 is None or C0 is None:        rng = np.random.default_rng(random_state)        A = rng.normal(size=(I, rank))        B = rng.normal(size=(J, rank))        C = rng.normal(size=(K, rank))    else:        A = _validate_factor(A0, I, rank, "A0")        B = _validate_factor(B0, J, rank, "B0")        C = _validate_factor(C0, K, rank, "C0")    B, _ = normalize_columns(B)    C, _ = normalize_columns(C)    weights = np.ones(rank, dtype=float)    history = []    prev = None    converged = False    n_iter = 0    for n_iter in range(1, int(max_iter) + 1):        A = update_A(X, B, C)        A, weights = normalize_columns(A)        B = update_B(X, A, C)        B, weights = normalize_columns(B)        C = update_C(X, A, B)        C, weights = normalize_columns(C)        Xhat = reconstruct(weights, A, B, C)        rec_error = relative_reconstruction_error(X, Xhat)        history.append(rec_error)        if prev is not None and abs(prev - rec_error) < float(tol):            converged = True            break        prev = rec_error    Xhat = reconstruct(weights, A, B, C)    return {        "weights": weights,        "A": A,        "B": B,        "C": C,        "Xhat": Xhat,        "E": X - Xhat,        "n_iter": n_iter,        "converged": converged,        "rec_error": relative_reconstruction_error(X, Xhat),        "history": history,    }  def parafac_als_multistart(X, rank, seeds, max_iter=100, tol=1e-8):    fits = [        parafac_als(X, rank, max_iter=max_iter, tol=tol, random_state=seed)        for seed in seeds    ]    best_index = int(np.argmin([fit["rec_error"] for fit in fits]))    return {        "fits": fits,        "best_index": best_index,        "best": fits[best_index],    } 

1. 1.

   Harshman, R. A. (1970). Foundations of the PARAFAC procedure: Models and conditions for an "explanatory" multi-modal factor analysis. UCLA Working Papers in Phonetics, 16, 1-84.

2. 2.

   Carroll, J. D., & Chang, J. J. (1970). Analysis of individual differences in multidimensional scaling via an N-way generalization of "Eckart-Young" decomposition. Psychometrika, 35(3), 283-319.

   doi:10.1007/BF02310791
3. 3.

   Kruskal, J. B. (1977). Three-way arrays: Rank and uniqueness of trilinear decompositions, with application to arithmetic complexity and statistics. Linear Algebra and its Applications, 18(2), 95-138.

   doi:10.1016/0024-3795(77)90069-6
4. 4.

   Bro, R. (1997). PARAFAC. Tutorial and applications. Chemometrics and Intelligent Laboratory Systems, 38(2), 149-171.

   doi:10.1016/S0169-7439(97)00032-4
5. 5.

   Andersson, C. A., & Bro, R. (2000). The N-way Toolbox for MATLAB. Chemometrics and Intelligent Laboratory Systems, 52(1), 1-4.

   doi:10.1016/S0169-7439(00)00071-X
6. 6.

   Sidiropoulos, N. D., & Bro, R. (2000). On the uniqueness of multilinear decomposition of N-way arrays. Journal of Chemometrics, 14(3), 229-239.

   doi:10.1002/1099-128X(200005/06)14:3<229::AID-CEM587>3.0.CO;2-N
7. 7.

   Bro, R., & Kiers, H. A. L. (2003). A new efficient method for determining the number of components in PARAFAC models. Journal of Chemometrics, 17(5), 274-286.

   doi:10.1002/cem.801
8. 8.

   Tomasi, G., & Bro, R. (2006). A comparison of algorithms for fitting the PARAFAC model. Computational Statistics & Data Analysis, 50(7), 1700-1734.

   doi:10.1016/j.csda.2004.11.013
9. 9.

   Stedmon, C. A., & Bro, R. (2008). Characterizing dissolved organic matter fluorescence with parallel factor analysis: a tutorial. Limnology and Oceanography: Methods, 6, 572-579.

   doi:10.4319/lom.2008.6.572
10. 10.

    Kolda, T. G., & Bader, B. W. (2009). Tensor Decompositions and Applications. SIAM Review, 51(3), 455-500.

    doi:10.1137/07070111X
11. 11.

    Kossaifi, J., Panagakis, Y., Anandkumar, A., & Pantic, M. (2019). TensorLy: Tensor Learning in Python. Journal of Machine Learning Research, 20(26), 1-6.