A Normal Form Algorithm for Tensor Rank Decomposition

We formulate our results for the case \(d = 3\), as the general case can be handled using the standard reshaping trick (Section 2.2). Our main contribution is a new numerical algorithm for computing the rank-\(r\) CPD of a third-order tensor \(\mathcal {A}\) based on linear algebra, under the assumption that \(\mathcal {A}\) is \(r\)-identifiable and a flattening of \(\mathcal {A}\) has rank \(r\), see assumption 2.1. We call this algorithm cpd_hnf. It is based on a well-known reformulation of the problem as a system of polynomial equations, which we solve using state-of-the-art methods from computational algebraic geometry. This results in 3.1. We show that this algorithm generalizes pencil-based algorithms [Beltrán et al. 2019; Leurgans et al. 1993; Lorber 1985; Sanchez and Kowalski 1990; De Lathauwer 2006] from very low ranks to much higher ranks in the unbalanced case; see theorem 3.13.23.3.

We give a new, explicit description of the complexity of the tensor rank decomposition problem in terms of an algebraic property of aforementioned polynomial system: the multigraded regularity of a non-saturated ideal in the homogeneous coordinate ring of \(\mathbb {P}^m \times \mathbb {P}^n\) (proposition 3.1). We characterize the regularity in terms of the rank of a structured matrix obtained directly from the rank-1 terms in Equation (CPD), see example 4.1. These new insights allow us to formulate a conjecture regarding the regularity (conjecture 1). We prove this conjecture for many formats, see theorem 4.2. These results are of independent interest for the field of polynomial system solving. They have the following consequence related to the complexity of our algorithm.

Our numerical experiments in Section 6 show that the proposed algorithm is highly efficient and gives accurate results (see Figure 2). For instance, we can compute the decomposition of a \(7 \times 7 \times 7 \times 7 \times 6 \times 6 \times 5 \times 5\) tensor of rank 1,000 in double-precision arithmetic with an accuracy of order \(10^{-15}\) in 441 s—a feat we believe has not been matched by other tensor decomposition algorithms. Moreover, cpd_hnf seems to behave well in the presence of noise.

A Julia implementation of cpd_hnf, including driver routines to reproduce our numerical experiments, is provided at https://gitlab.kuleuven.be/u0072863/homogeneous-normal-form-cpd.

The idea of computing tensor decompositions via polynomial root finding is central in apolarity-based approaches such as Bernardi et al. [2011];, 2013]; Bernardi and Taufer [2020], and Brachat et al. [2010]; Nie [2017] for the symmetric case. The Hankel operators play the role of normal forms in this context. These operators can be obtained partially from the tensor \(\mathcal {A}\). In most cases, an additional polynomial system needs to be solved to complete the Hankel operators [Mourrain 2018, Section 4]. For instance, this is step (2) in Algorithm 5.1 of Brachat et al. [2010]. Although our method works only under certain assumptions on the rank of \(\mathcal {A}\), in contrast to apolarity-based methods it requires only linear algebra computations, and it operates in polynomial rings with fewer variables. Moreover, the method from Brachat et al. [2010] uses connected-to-one bases for its normal form computations. The choice of such a basis is discussed at length in Sections 4.3 and 7.2 of Bernardi and Taufer [2020]. In this article, we exploit the flexibility of truncated normal forms [Telen et al. 2018] to achieve better numerical results. For a comparison, see Section 4.3.3 of Telen [2020b].

Kuo and Lee [2018] obtained an affine system of polynomial equations as in example 1.1, but it is solved using homotopy continuation methods. This approach is infeasible for some formats that are handled without any problems by our algorithm. For instance, for the aforementioned eighth-order tensor of rank 1,000, the method of Kuo and Lee [2018] would need to track over \(40 \times 10^9\) homotopy paths. We argue that the eigenvalue-based methods proposed in this article are more natural to use in an overdetermined setting.

The state-of-the-art algorithms for tensor rank decomposition using only linear algebra computations were proposed by Domanov and De Lathauwer [2014];, 2017]. Although these methods work under slightly milder conditions, our numerical experiments suggest that Algorithm 1 of Domanov and De Lathauwer [2017] often requires the construction of larger matrices than those in our algorithm. There is no explicit connection with polynomial equations. The algorithm and its complexity depend on a parameter \(l\), which is chosen incrementally by trial and error for each format. In an analogous way, the complexity of our algorithm is governed by the choice of a parameter. However, our analysis in Section 4 tells us a priori which parameter value should be used, circumventing a trial-and-error approach. In Section 6, we demonstrate that cpd_hnf improves on Domanov and De Lathauwer [2014];, 2017] in terms of computational complexity and accuracy.

In Section 2, we state our assumptions and show how computing the CPD of a tensor \(\mathcal {A}\) is formulated as a system of polynomial equations. In Section 3, we make the connection with normal form methods explicit. That is, we describe how a pre-normal form can be computed directly from the tensor \(\mathcal {A}\) and how this allows us to reduce the above polynomial system to an eigenvalue problem. We explain how the approach generalizes so-called pencil-based algorithms in Section 3.3. This leads to a complete, high-level description of our algorithm cpd_hnf in Section 3.4. In Section 4, we study the regularity of the ideal associated to our polynomial system. These results are the starting point for our analysis of the complexity of cpd_hnf in Section 5. In Section 6, we demonstrate the efficiency and accuracy of cpd_hnf relative to the state of the art through several numerical experiments. The article is completed in Section 7 with some final conclusions.

For the moment, let us assume \(\mathcal {A} \in \mathbb {C}^{\ell +1}\otimes \mathbb {C}^{m+1}\otimes \mathbb {C}^{n+1}\) is a third-order tensor of rank \(r\). The cpd_hnf algorithm works under the following assumption.

By generic, we mean that \(\mathcal {A}\) is contained in a Zariski dense open subset of the set of all rank-\(r\) tensors. We will describe this open subset more explicitly in lemma 2.2. For now, we point out that in this open set, each tensor \(\mathcal {A}\) is such that

These are necessary conditions for our algorithm to work. We now briefly recall their meaning.

Condition (i) is usually very weak and only of a technical nature. Recall that the set of all rank-1 tensors forms an algebraic variety, i.e., the solution set of a system of polynomial equations, called the Segre variety \(\mathcal {S}\). A rank-\(r\) CPD of a tensor \(\mathcal {A}\) is a set of \(r\) rank-1 tensors whose sum is \(\mathcal {A}\). The set of all such rank-\(r\) CPDs is denoted by \(\mathcal {S}^{[r]} = \lbrace \mathcal {X} \subset \mathcal {S} \mid |\mathcal {X}|=r \rbrace\). Tensor rank decomposition consists of computing an element of the fiber ofat a rank-\(r\) tensor \(\mathcal {A}\). For brevity, we denote the image of \(f\) byTensor \(\mathcal {A}\) is called \(r\)-identifiable if the fiber \(f^{-1}(\mathcal {A})\) contains exactly one element. That is, there is a unique set in \(\mathcal {S}^{[r]}\) whose sum is \(\mathcal {A}\). Generally, \(r\)-identifiability fails only on a strict closed subvariety of (the Zariski closure of) \(\mathcal {S}_r\); see Chiantini and Ottaviani [2012]; Chiantini et al. [2014];, 2017]; Bocci et al. [2014]. This property is called the generic \(r\)-identifiability of \(\mathcal {S}\). These results entail that \(r\)-identifiability fails only on a subset of \(\mathcal {S}_r\) of Lebesgue measure zero if the rank \(r\) and dimensions \((\ell +1,m+1,n+1)\) satisfy some very weak conditions; see Section 1 in Chiantini et al. [2014] and Section 3 in Chiantini et al. [2017] for more details and statements for higher-order tensors as well. If \(\mathcal {S}_r\) is generically \(r\)-identifiable, then there is a Zariski-open subset of the closure of \(\mathcal {S}_r\) so that \(f^{-1}\) is an analytic, bijective tensor decomposition function. We aspire to solve the tensor decomposition problem only in this well-behaved setting.

Condition (ii) is more restrictive, but allows us to tackle the tensor decomposition problem using only efficient linear algebra. Recall that the standard 1-flattening of \(\mathcal {A} \in \mathbb {C}^{\ell +1} \otimes \mathbb {C}^{m+1} \otimes \mathbb {C}^{n+1}\) consists of interpreting \(\mathcal {A}\) as the matrix \(\mathcal {A}_{(1)} \in \mathbb {C}^{(\ell +1) \times (m+1)(n+1)}\). For a rank-1 tensor \(\mathcal {A}=\alpha \otimes \beta \otimes \gamma\) this identification is defined byand the general case follows by linearity. The tensor product² in the foregoing expression is also called the reverse-order Kronecker productwhere the standard bases were assumed for these Euclidean spaces and the indices \((i,j)\) are sorted by the reverse lexicographic order. Note that the 1-flattening is easy to compute when \(\mathcal {A}\) is given in coordinates relative to the standard tensor product basis. In that case it suffices to reshape the coordinate array to an \((\ell +1) \times (m+1)(n+1)\) array (e.g., as in Julia’s or Matlab’s reshape function).

For dealing with higher-order tensors \(\mathcal {A}\in \mathbb {C}^{n_1+1}\otimes \cdots \otimes \mathbb {C}^{n_d+1}\) with \(d \gt 3\), we rely on the well-known reshaping trick. It consists of interpreting a higher-order tensor as a third-order tensor. While the approach described in the remainder of the article could be applied directly to higher-order tensors as well, the range of ranks \(r\) to which this version would apply is (much) more restrictive than the range obtained from reshaping.

Recall that reshaping \(\mathcal {A}\) to a third-order tensor \(\mathcal {A}_{(I,J,K)}\) is a linear isomorphismthat identifies \(\alpha ^1 \otimes \cdots \otimes \alpha ^d\) with \((\otimes _{i\in I} \alpha ^i) \otimes (\otimes _{j\in J} \alpha ^j) \otimes (\otimes _{k\in K} \alpha ^k)\), wherein \(I \sqcup J \sqcup K\) partitions \(\lbrace 1,\ldots ,d\rbrace\). The general case follows by linearity from the universal property [Greub 1978].

It was shown in Section 7 of Chiantini et al. [2017] that under some conditions, the unique CPD of \(\mathcal {A}\) can be recovered from the CPD of a reshaping \(\mathcal {A}_{(I,J,K)}\). Here, we exploit the following slightly more general result that requires no \(r\)-identifiability of \(\mathcal {A}\).

For decomposing higher-order tensors \(\mathcal {A}\) of rank \(r\) via the reshaping trick, we proceed as follows:

Note that the rank-1 decompositions in the third step can be computed for example with a sequentially truncated higher-order singular value decomposition [Vannieuwenhoven et al. 2012], as in our implementation, or with a cross approximation [Oseledets et al. 2008].

Because of the reshaping trick, we will henceforth describe our approach only for \(d=3\).

Having delineated the range (R) in which the proposed cpd_hnf algorithm will work, we continue by describing it. Our strategy to compute the tensor rank decomposition ofis to compute the pointsfirst as the solution of a system of (low-degree) polynomial equations. Thereafter, we compute the \(\alpha _{i}\)’s with efficient linear algebra by plugging \(\beta _{i}\) and \(\gamma _{i}\) into Equation (A) and solving the resulting linear system of equations. Indeed, sinceand assumption 2.1 guarantees that the second \(r \times (m+1)(n+1)\) matrix is of full rank \(r\) (see lemma 2.2 below), so it has a right inverse. Applying the latter on the right to \(\mathcal {A}_{(1)}\) yields the corresponding \(\alpha _i\)’s. Note that it suffices to compute the points \(\beta _{i}\) and \(\gamma _{i}\) up to a nonzero scaling factor. Therefore, it is natural to consider our problem in complex projective space.

Recall that the \(k\)-dimensional complex projective space \(\mathbb {P}^{k}\) is the space of equivalence classes \([x] = \lbrace \lambda x \mid \lambda \in \mathbb {C} \setminus \lbrace 0\rbrace \rbrace\) for \(x \in \mathbb {C}^{k+1} \setminus \lbrace 0\rbrace\). The entries of the vector \(x = (x_0, \ldots , x_k) \in \mathbb {C}^{k+1} \setminus \lbrace 0\rbrace\) are called homogeneous coordinates of \([x]\). With a slight abuse of notation, we will write \(x = (x_0 : \cdots : x_k) \in \mathbb {P}^k\) for both the equivalence class \([x]\) and a set of homogeneous coordinates \(x \in \mathbb {C}^{k+1} \setminus \lbrace 0\rbrace\).

The proposed cpd_hnf algorithm exploits the fact that the points \(\zeta _i = (\beta _i, \gamma _i) \in \mathbb {P}^m \times \mathbb {P}^n\) are defined by algebraic relations on \(X = \mathbb {P}^{m} \times \mathbb {P}^n\) that can be computed directly from the tensor \(\mathcal {A}\). Such algebraic relations are homogeneous polynomials in a bi-graded polynomial ring. In this context, it is natural to think of degrees as 2-tuples \((d,e) \in \mathbb {N}^2\), where \(d\) is the degree in the variables corresponding to \(\mathbb {P}^m\) and \(e\) is the degree in the variables corresponding to \(\mathbb {P}^n\). Note that this differs from the more familiar setting where the degree is a natural number \(d \in \mathbb {N}\). For more on multi-graded rings, see Miller and Sturmfels [2005]. The flexibility provided by this bi-graded setting reduces the complexity of our algorithm. Concretely, we work in the \(\mathbb {N}^2\)-graded polynomial ringHere, we used the notation \(x^{a} = x_{0}^{a_{0}} \cdots x_{m}^{a_{m}}\) and \(|a| = a_{0} + \cdots + a_{m}\) for \(a = (a_{0}, \ldots , a_{m}) \in \mathbb {N}^{m + 1}\), and analogously for \(b\). The graded pieces \(S_{(d,e)}\) are vector spaces over \(\mathbb {C}\). The variables \(x\) correspond to homogeneous coordinates on the factor \(\mathbb {P}^{m}\) and in \(S_{(d,e)}\), \(d\) is the degree in the \(x\)-variables. Analogously, elements of \(S_{(d,e)}\) have degree \(e\) in the \(y\)-variables, which are related to \(\mathbb {P}^n\).

An element \(f \in S\) is called homogeneous if \(f \in S_{(d,e)}\) for some \((d,e)\in \mathbb {N}^2\). The key reason that the ring \(S\), with its grading (2.2), is naturally associated to \(X\) is that homogeneous elements \(f \in S_{(d,e)}\) have well-defined zero sets on \(X\). By this, we mean that for \(f \in S_{(d,e)}\) and for any \(\zeta = (x, y) \in X\), \(f(x,y) = 0\) is independent of the choice of homogeneous coordinates. Indeed, this follows fromTherefore, whenever \(f\) is homogeneous it makes sense to write \(f(\zeta) = 0\) if \(f\) vanishes on some set of homogeneous coordinates for \(\zeta\), and to define the subvariety of \(X\) defined by \(f\) as \(V_X(f) = \lbrace \zeta \in X ~|~ f(\zeta) = 0\rbrace\). For a homogeneous ideal \(I \subset S\) (i.e., \(I\) is generated by homogeneous elements with respect to the grading (2.2)), we denote the subvariety of \(X\) corresponding to \(I\) byThat is, \(V_X(I)\) contains the solutions of the polynomial system defined by the simultaneous vanishing of all homogeneous equations \(f \in I\).

With this notation in place, we turn back to the tensor \(\mathcal {A}\) and show that, under suitable assumptions, \(V_X(I) = \lbrace \zeta _1,\ldots ,\zeta _r\rbrace\) for some homogeneous ideal \(I\) defined by \(\mathcal {A}\). This will generalize the procedure in example 1.1. Consider the flattening \(\mathcal {A}_{(1)}\) from Equation (2.1). This flattening has a natural interpretation as a \(\mathbb {C}\)-linear mapIndeed, we can identify the space \(\mathbb {C}^{m+1} \otimes \mathbb {C}^{n+1}\) with the graded piece \(S_{(1,1)}\) of degree \((1,1)\) bywhere \(e_{k}\) represents the \((k+1)\)-st standard basis vector of \(\mathbb {C}^{m + 1}\), and analogously for \(e_{l}\) and \(\mathbb {C}^{n+1}\). For brevity, we also write \(\mathcal {A}_{(1)}\) for a matrix representation of the map \(\mathcal {A}_{(1)}\), where the standard basis of monomials (2.4) is used for \(S_{(1,1)}\). The \(\ell +1\) rows of \(\mathcal {A}_{(1)}\) are elements of the dual space \((S_{(1,1)})^\vee\) represented in its dual basisFrom this vantage point, the rows of the second factor in Equation (2.1) areIt is clear that for any \(f \in S_{(1,1)}\), we have \((v \otimes w)^\top (f) = f(v,w)\).

Let \(f_1, \ldots , f_s \in S_{(1,1)}\) be a \(\mathbb {C}\)-basis for the kernel \(\ker \mathcal {A}_{(1)}\). The \(f_i\) generate a homogeneous ideal \(I = \langle {\ker \mathcal {A}_{(1)}} \rangle = \langle {f_1,\ldots ,f_s} \rangle \subset S\). The following lemma characterizes precisely what we mean by generic in assumption 2.1.

Point (iii) in lemma 2.2 is a technicality we will need in Section 3 to invoke the results from Section 5.5 in Telen [2020b]. For the reader who is familiar with algebraic geometry, we included a consequence of lemma 2.2 relating \(I\) to the vanishing ideal of \(V_X(I)\) in corollary A.1.

We were not able to construct an example for which (i)+(ii) are satisfied, but (iii) is not. This raises the question whether (i)+(ii) implies (iii), and if not, how our method extends to counterexamples.

We conclude the section by presenting a high-level pseudo-algorithm for computing the tensor rank decomposition (A) of \(\mathcal {A}\). This is presented as 2.1. Its steps summarize the discussion up to this point.

The main focus of this article is dealing with step 4 of 2.1. We will employ a state-of-the-art method for solving polynomial equations on \(X\), based on homogeneous normal forms. This strategy is described in the next section.

Our algorithm is built on a multi-homogeneous version of the eigenvalue theorem (theorem 3.1 below), which allows us to find solutions of systems of polynomial equations via eigenvalue computations. Behind this is the theory of homogeneous normal forms. In our context, these are \(\mathbb {N}^2\)-graded versions of truncated normal forms, introduced in Telen et al. [2018], and special cases of the more general toric homogeneous normal forms used by Telen [2020a]; Bender and Telen [2020] and formally introduced in Telen [2020b]. For our purpose, it suffices to work with slightly simpler objects, called pre-normal forms, so we use homogeneous normal forms only implicitly. For full proofs and more details, we refer to Section 5.5.4 in Telen [2020b].

Consider the space \(X = \mathbb {P}^m \times \mathbb {P}^n\) and its associated ring \(S\), which is the \(\mathbb {N}^2\)-graded polynomial ring from Equation (2.2). Let \(I \subset S\) be a homogeneous ideal. The ideal \(I\) inherits the grading on \(S\):is a \(\mathbb {C}\)-subspace of \(S_{(d,e)}\), and the quotient ring(these are quotients of vector spaces) similarly inherits this grading. Our objective is to compute a set of homogeneous coordinates of the points in \(V_X(I)\) by using linear algebra computations. To accomplish this, it is necessary to work in graded pieces of \(S\), \(I\), and \(S/I\). Let \(M\) denote either of the latter. The (multi-graded) Hilbert function \(\mbox{HF}_M : \mathbb {N}^2 \rightarrow \mathbb {N}\) is given byand keeps track of the dimension of the vector space \(M_{(d,e)}\). Note that for a homogeneous ideal \(I \subset S\) and \((d,e) \in \mathbb {N}^2\), we have \(\mbox{HF}_{S/I}(d,e) = \mbox{HF}_S(d,e) - \mbox{HF}_I(d,e)\).

Before stating the main result of this section, theorem 3.1, we present three auxiliary lemmas. Lemma 3.13.2 are well-known in commutative algebra. For lack of a precise reference and completeness, we included short proofs in the Appendix.

Linear maps \(N: S_{(d,e)} \rightarrow \mathbb {C}^r\) satisfying Equation (3.1) will play an important role in translating our polynomial root finding problem into an eigenvalue problem. We will therefore give them a name.

If \(r = \ell +1\), then assumption 2.1 implies that \(\mathcal {A}_{(1)}: S_{(1,1)} \rightarrow \mathbb {C}^r\) is a pre-normal form on \(S_{(1,1)}\).

We use the terminology pre-normal form, because, following Definition 5.5.6 in Telen [2020b], the term normal form is reserved for linear maps satisfying an extra technical condition. Normal forms are closely related to a classical result in computational algebraic geometry called the eigenvalue theorem [Cox et al. 2006, Section 2.4], which is used for computing isolated solutions to polynomial systems. The pre-normal forms introduced above will be useful for formulating a multi-homogeneous version of this theorem, namely, theorem 3.1 below. Before stating it, we need to fix some additional notation.

In lemma 3.3, \(N_{h_0}\) is the composition of \(N\) with the linear map that represents multiplication by \(h_0 \in S_{(d^{\prime },e^{\prime })}\). More generally, for \(g \in S_{(d^{\prime },e^{\prime })}\), we define \(N_g : S_{(1,1)} \rightarrow S_{(d,e)}\) as \(N_g(f) = N(gf)\). In theorem 3.1, we will restrict the map \(N_{h_0}\) to an \(r\)-dimensional subspace \(B \subset S_{(1,1)}\) such that the resulting map is invertible. We denote this restriction by \((N_{h_0})_{|B}\). In practice, it suffices to select \(r\) columns of the matrix of \(N_{h_0}\) so that the resulting \(r \times r\) submatrix \((N_{h_0})_{|B}\) is invertible.

In what follows, we will write \(w_i ^\top = (\beta _i \otimes \gamma _i)^\top _{|B} \in B^\vee\) for the linear functional representing evaluation at \(\zeta _i = (\beta _i, \gamma _i)\). This is the restriction of the functional \((\beta _i \otimes \gamma _i)^\top \in S_{(1,1)}^\vee\) to the vector space \(B\). Concretely, we setAfter fixing a basis for \(B^\vee\), \(w_i^\top\) can be represented as a row vector.

If \(g, h \in S_{(d^{\prime },e^{\prime })} \setminus \lbrace 0\rbrace\) are two homogeneous polynomials of the same degree, then the fraction \(g/h\) is a well-defined function on \(X \setminus V_X(h)\). Indeed, the evaluation of this function does not depend on the choice of homogeneous coordinates. Therefore, we may write \((g/h)(\zeta)\) for the evaluation of this rational function at \(\zeta \in X \setminus V_X(h)\).

After fixing a basis for \(B\) and representing \(w_i\) in the dual basis for \(B^\vee\), Equation (3.3) is a standard matrix eigenproblem: \(w_i^\top M_{g/h_0} = \lambda _i w_i^\top\). That is, \((\lambda _i,w_i)\) is a left eigenpair of the \(r\) \(\times\) \(r\) matrix \(M_{g/h_0}\). Note that theorem 3.1 implies that all maps of the form \(M_{g/h_0}\) share a set of eigenvectors.

We now sketch one way of using theorem 3.1 to retrieve the coordinates of \(\zeta _i = (\beta _i, \gamma _i)\) from eigenvalues, assuming a pre-normal form \(N: S_{(d,e)} \rightarrow \mathbb {C}^r\) is given. The problem of computing a pre-normal form is addressed in the next subsection. We assume \(d \ge 2\).⁴ Let \(h \in S_{(d^{\prime }-1,e^{\prime })}\) and \(h_0 \in S_{(d^{\prime },e^{\prime })}\) be homogeneous polynomials that do not vanish at any of the points \(\zeta _i\). These can be chosen generically. Set \(g_j = x_j h \in S_{(d^{\prime },e^{\prime })}, j = 0, \ldots , m\). Choose \(B \subset S_{(1,1)}\) of dimension \(r\) such that \((N_{h_0})_{|B}\) is invertible and compute the matrices \(M_j = M_{g_j/h_0} = (N_{h_0})_{|B}^{-1} \circ (N_{g_j})_{|B}\). By theorem 3.1, the eigenvalues of \(M_j\) are given by \(\lambda _{ji} = (g_j/h_0)(\zeta _i)\). Writing \(\beta _{ij}\) for the \(j\)th coordinate of \(\beta _i\), we haveNote that if \((d,e) = (2,1)\), we can take \(h = 1\).

Subsequently, we compute \(\gamma _i\) by solving the linear system of equations

The foregoing approach requires that \((d,e) \ne (1,1)\). Otherwise \(S_{(d^{\prime },e^{\prime })}=S_{(0,0)} = \mathbb {C}\), and we can only evaluate constant functions using theorem 3.1.

As illustrated in the previous subsection, once we have computed a pre-normal form, the points \(\zeta _i\) can be recovered using basic linear algebra computations. A natural next issue to address is how to compute a pre-normal form.

Our starting point is a basis \(f_1, \ldots , f_s \in S_{(1,1)}\) of \(\ker \mathcal {A}_{(1)}\), generating our ideal \(I = \langle {\ker \mathcal {A}_{(1)}}\rangle = \left\langle {f_1,\ldots ,f_s} \right\rangle\). For any tuple \((d,e) \in \mathbb {N}^2\) such that \((d,e) \ge (1,1)\), the degree-\((d,e)\) part \(I_{(d,e)}\) of \(I\) is the \(\mathbb {C}\)-vector space spanned bywhere \((d^{\prime },e^{\prime }) = (d-1,e-1)\). If \(d = 0\) or \(e = 0\), then we have \(I_{(d,e)} = \lbrace 0\rbrace\). In analogy with Equation (1.1), we construct a matrix whose rows are indexed by the monomial basis elements of \(S_{(d,e)}\) (i.e., the monomials \(\lbrace x^ay^b ~|~ |a| = d, |b|=e \rbrace\)), and whose columns are the polynomials (3.5) expanded in this basis. We denote this matrix bySuch matrices represent graded resultant maps in the terminology of Section 5.5.4 in Telen [2020b]. They are multihomogeneous versions of the classical Macaulay matrices [Macaulay 1916]. We present an explicit example below in example 3.2.

Observe that the Hilbert function \(\mbox{HF}_I(d,e)\) is given by the rank of \(R_I(d,e)\) and \(\mbox{HF}_{S/I}(d,e)\) by its corank. This follows immediately from the observation that the columns of \(R_I(d,e)\) span \(I_{(d,e)}\). A left nullspace matrix \(N\) of \(R_I(d,e)\) represents a map \(S_{(d,e)} \longrightarrow S_{(d,e)}/I_{(d,e)} \simeq \mathbb {C}^{\mbox{HF}_{S/I}(d,e)}\). This has the following consequence.

We conclude that a pre-normal form \(N\) can be computed, for instance, from a full singular value decomposition (SVD) of \(R_I(d,e)\), where \(\mbox{HF}_{S/I}(d,e) = r\). This solves the problem of computing a pre-normal form, assuming that we know a degree \((d,e) \in \mathbb {N}^2\) for which \(\mbox{HF}_{S/I}(d,e) = r\). The problem of finding such degrees is addressed in Section 4.

To obtain the homogeneous coordinates for \(\zeta _1, \ldots , \zeta _r\) as eigenvalues of the matrices \(M_{g/h_0}\), we usually have to work with pre-normal forms on \(S_{(d,e)}\), where \((d,e) \ne (1,1)\) and \((d^{\prime },e^{\prime }) \ge (0,0)\). An exception is the case where \(r \le m+1 \le \ell + 1\). For these tensors of very low rank, a pre-normal form \(N: S_{(1,1)} \rightarrow \mathbb {C}^r\) will suffice under the mild condition thatThe underlying reason is that vanishing at \(\lbrace \zeta _1,\ldots , \zeta _r\rbrace\) gives \(r\) linearly independent conditions on \(S_{(1,0)}\). The proof of the following theorem is another consequence of the theory of homogeneous normal forms and is deferred to the Appendix.

This theorem is exploited to compute \(\zeta _i = (\beta _i, \gamma _i)\) efficiently as follows. If \(\mathcal {A}\) satisfies assumption 2.1 and Equation (3.7), then we take a basis of the \(r\)-dimensional row span of \(\mathcal {A}_{(1)}\) in theorem 3.2 as our pre-normal form. This can be obtained from a compact SVD of \(\mathcal {A}_{(1)}\). Once the \(\gamma _i\) are computed from the eigenvalues of the \(M_{y_j/h_0}\), the \(\beta _i\) can be obtained as in Equation (3.4). Alternatively, one can use the eigenvectors of these commuting matrices \(M_{y_j/h_0}\) for \(j=1,\ldots ,r\) [Telen 2020b, Theorem 5.5.3].

Theorem 3.2 is intimately related to what Beltrán et al. [2019] called pencil-based algorithms for solving Equation (CPD) when the rank satisfies \(r \le m+1 \le \ell +1\), such as those by Leurgans et al. [1993]; Lorber [1985]; Sanchez and Kowalski [1990]. Recall that pencil-based algorithms assume that \(\mathcal {A}\in \mathbb {C}^{r \times r \times (n+1)}\) is a rank-\(r\) tensor.⁵ In addition, they assume that the \(\alpha _i\) form a linearly independent set, and likewise for the \(\beta _i\)’s. Then, we have that the tensor contraction of \(\mathcal {A}\) with \(h_0^\top \in (\mathbb {C}^{n+1})^\vee\), i.e.,is an invertible \(r \times r\) matrix insofar as \(h_0^\top \gamma _i \ne 0\). Herein, \(A \in \mathbb {C}^{r \times r}\) (respectively, \(B \in \mathbb {C}^{r \times r}\)) has the \(\alpha _i\)’s (respectively, \(\beta _i\)’s) as columns, and \(D_{h_0} = \mbox{diag}(h_0^\top \gamma _1, \ldots , h_0^\top \gamma _r)\). Let \(\widetilde{N_{h_0}} = h_0^\top \cdot _3 \mathcal {A}\) and \(\widetilde{N_{g}} = g^\top \cdot _3 \mathcal {A}\) for \(h_0, g \in (\mathbb {C}^{n+1})^\vee\). Then, we haveso that the points \(\beta _i\) can be recovered uniquely from the matrix of eigenvectors \(B^\top\), provided that \(h_0^\top \gamma _i \ne 0\) for all \(i=1,\ldots ,r\). The \(\alpha _i\)’s and \(\gamma _i\)’s can then be recovered from the 2-flattening; see Leurgans et al. [1993]; Lorber [1985]; Sanchez and Kowalski [1990]; Beltrán et al. [2019]; De Lathauwer [2006] for more details. With the foregoing suggestive notation, it is easy to see that the matrix of \(\widetilde{N_{h_0}} : f \mapsto N(fh_0)\) with respect to the standard bases is precisely Equation (3.8). Indeed, note that since we can take \(N = \mathcal {A}_{(1)}\), we have \(f h_0 \simeq f \otimes h_0\) and so \(N(f h_0) = \mathcal {A}_{(1)}(f \otimes h_0)\).

Pencil-based algorithms may thus be interpreted as a special case of the proposed cpd_hnf algorithm based on homogeneous normal forms when \(r \le m+1 \le \ell +1\). Note that because of the numerical instabilities analyzed by Beltrán et al. [2019] caused by extracting \(\beta _i\) from the eigenvectors, we prefer to extract the \(\zeta _i = (\beta _i,\gamma _i)\) in a different way. We compute \(\beta _i\) from the eigenvalues of \(M_{x_i/h_0}\) and the corresponding \(\gamma _i\) from the linear system Equation (3.4).

The discussion so far is distilled into 3.1. This algorithm implements step 4 of 2.1. Note that we dropped the tilde on top of the \(N_\star\)’s in lines 2–7 to streamline the presentation.

The first phase of the algorithm, up to line 21, constructs the pre-normal form \(N\) and chooses an \(N_{h_0}\). This phase depends on whether we can invoke the more efficient theorem 3.2 (\(r \le m+1\)) or we need the full power of theorem 3.1. In the former case, we can take \(N = \mathcal {A}_{(1)}\), while in the latter case, we need to take \(N\) equal to the left null space of \(R_I(d,e)\). How we choose the degree \((d,e)\) in line 9 is explained in Section 4. The matrix \(R_I(d,e) \in \mathbb {C}^{\mbox{HF}_S(d,e) \times s \mbox{HF}_S(d^{\prime },e^{\prime })}\) can be constructed efficiently column-by-column without polynomial multiplication. Indeed, by Equation (3.5) it suffices to copy the coefficients of \(f_i\) relative to the monomial basis of \(S_{(1,1)}\) into the correct rows; see also example 3.2 below. The left null space \(N\) can be extracted from the last \(r\) columns of the \(U\)-factor in the SVD \(R_I(d,e)=USV^H\), where \(\cdot ^H\) denotes the conjugate transpose. In our implementation, the matrix \(N_{h_0} \in \mathbb {C}^{r \times \mbox{HF}_S(1,1)}\) is chosen by sampling the coefficients of \(h_0 \in S_{(d^{\prime },e^{\prime })}\) identically and independently distributed from a Gaussian distribution. With probability 1, \(h_0\) satisfies \(h_0(\gamma _i) \ne 0\) for all \(i\); hence, this is a valid choice of \(h_0\).

The next phase of the algorithm, in lines 21–22, chooses a basis \(B\). Although in theory theorem 3.1 enables us to choose any \(B\) such that \((N_{h_0})_{|B}^{-1}\) is invertible, Telen et al. [2018] showed that for reasons of numerical stability it is crucial to choose \(B\) such that \((N_{h_0})_{|B}\) is well-conditioned. In practice, such a subspace \(B\) can be found using a QR decomposition with optimal column pivoting or by using the SVD [Telen 2020b, Chapter 4]. We stated the QR approach in 3.1.

The multiplication matrices are constructed straightforwardly as the formula suggests in lines 23 to 25. Note that the upper triangular matrix \((N_{h_0})_{|B}\) does not need to be inverted explicitly, rather the system can be solved by backsubstitution.

In line 26, the matrices \(M_{(h x_k)/h_0}\) need to be simultaneously diagonalized, as we have that \(M_{(h x_k)/h_0} = V^{-1} \operatorname{diag}(\beta _k) V\). We compute \(V\) from a random linear combination of the \(M_{(h x_k)/h_0}\)’s, and then extract \(\beta _k\) as the diagonal entries from the (approximately) diagonalized matrix \(V^{-1} M_{(h x_k)/h_0} V\).

The system in line 28 is solved efficiently by noting that the coefficients of \(f_j\) can be arranged in a matrix \(F_j\) of size \((m+1) \times (n+1)\), so that \(f_j(x,y) = x^\top F_j y\). Hence, Equation (3.4) boils down to computing the kernel of \(Ay = 0\) where the rows of \(A \in \mathbb {C}^{s \times (n+1)}\) are the row vectors \(\beta _i^\top F_j\). Note that by assumption 2.1, \(\ker A\) is spanned by \(\gamma _i\). In line 29 the obtained solution \((\beta _i,\gamma _i)\) is refined using standard Newton iterations. Here the pseudo-inverse of the Jacobian matrix of \(f_1, \ldots , f_s\) is used. This matrix has full rank if and only if \((\beta _i, \gamma _i)\) has multiplicity one. This can be used to check condition (iii) in lemma 2.2.

We now present two illustrative examples. The first one shows how to use the techniques explained above on the tensor \(\mathcal {A}\) in example 1.1.

The ideal \(I\) in the previous example has the property that \(\mbox{HF}_{S/I}(1+d^{\prime },1+e^{\prime }) = r\) for all \((d^{\prime },e^{\prime }) \in \mathbb {N}^2\). Our next example shows that this is not the case in general.

The algorithm is implemented for real and complex input tensors. In the former case, all computations are performed over the reals with the exception of the computation of the left null space of \(R_I(d,e)\). In the real case, the algorithm continues with the real part of the output of this step.

At the start of the algorithm, we compress the \((\ell +1) \times (m+1) \times (n+1)\) input tensor, which is to be decomposed into \(r\) rank-1 terms, to a \(\min \lbrace \ell +1,r\rbrace \times \min \lbrace m+1,r\rbrace \times \min \lbrace n+1,r\rbrace\) tensor. For this, we apply ST-HOSVD compression [Vannieuwenhoven et al. 2012] with truncation rank \((\min \lbrace \ell +1,r\rbrace , \min \lbrace m+1,r\rbrace , \min \lbrace n+1,r\rbrace)\). Most of the experiments below are chosen so that this step performs no computations.

The kernel of \(\mathcal {A}_{(1)}\) in 2.1 is computed by an SVD. The linear system at the end of 2.1 is solved by computing the Khatri–Rao product \(K = [\beta _i \otimes \gamma _i]_{i=1}^r\) and then solving the overdetermined linear system \(K A = \mathcal {A}_{(1)}^\top\) for \(A\). The rows of \(A\) then correspond to the \(\alpha _i\)’s.

In 3.1, the degree \((d,e)\) is determined automatically by assuming conjecture 1 is true and selecting a valid \((d,1)\) or \((1,e)\) based on a heuristic that takes into account the estimated computational cost and numerical considerations.¹²

The key computational bottleneck is the computation of the left nullspace of \(R_I(d,e)\) in line 11 of 3.1. When the number of entries of \(R_I(d,e)\) is smaller than 10,000, we use the standard SVD-based kernel computation. In larger cases, for efficiency, we propose to employ Arpack.jl’s eigs function to extract the left null space from the Hermitian Gram matrix \(G = R_I(d,e) (R_I(d,e))^H\), where \(^H\) denotes the Hermitian transpose. The eigs function (with parameters tol = 1e-6 and maxiter = 25) extracts the \(r\) eigenvalues of smallest modulus and their corresponding eigenvectors. Note that Arpack.jl employs a Bunch–Kaufman factorization [Bunch and Kaufman 1977] of the positive semidefinite input matrix \(G\) to perform its Lanczos iterations. Using eigs was up to \(75\%\) faster than the SVD-based approach for large problems.

The final step of 3.1 is implemented as discussed at the end of Section 3.4. The kernel is computed with an SVD and \(\gamma _i\) is taken as the left singular vector corresponding to the smallest singular value. After this step, we find \((\beta _i, \gamma _i)\). Observe that \(f_j(\beta _i, \gamma _i)\) should vanish exactly for all \(j = 1, \ldots , s\). We propose to refine the approximate roots \((\beta _i, \gamma _i)\) by applying three Newton iterations on \(f_1 = \cdots = f_s = 0\). This adds a minor computational overhead of \(r(m+1)(n+1)(m+n)^2\).

We compare and evaluate the impact on numerical accuracy of four variants of 3.1. Two options are considered in each combination:

We append _eigs, _eigs+newton, _svd, and _svd+newton to cpd_hnf to distinguish the variants.

The experimental setup is as follows. Random rank-\(r\) tensors are generated as explained at the start of this section, sampling one rank-\(r\) tensor for each combination \((\ell +1,m+1,n+1,r)\) with \(25 \ge m + 1 \ge n+1 \ge 2\) and \(\ell +1=r=\min \lbrace \mathcal {R}(m,n,(2,1)),mn\rbrace\) where \(\mathcal {R}\) is as in Equation (4.9). For each variant of cpd_hnf, the relative backward error\(\Vert \mathcal {A}- \sum _{i=1}^r \alpha _i^{\prime } \otimes \beta _i^{\prime } \otimes \gamma _i^{\prime } \Vert _F \Vert \mathcal {A}\Vert _F^{-1},\) and execution time is recorded, where \((\alpha _i^{\prime }, \beta _i^{\prime }, \gamma _i^{\prime })\) are the estimates obtained by the algorithm.

The experiment was performed on cascadelake nodes and the results are displayed in Figure 1. It is visually evident in the right panels that the Newton refinement improves the overall relative backward error for both eigs and svd by several orders of magnitude and brings both to seemingly the same level of about \(10^{-16}\)–\(10^{-12}\). However, there is a marked difference in accuracy between eigs and svd before Newton refinement, the svd being more accurate by about 1–3 orders of magnitude. It is interesting to compare these results with the leftmost panel of Figure 2(b), which shows the corresponding results on the same tensors for a conceptually similar state-of-the-art method from the literature. The method also requires a kernel computation that, in the displayed graph, is performed with a singular value decomposition. This means that the methods in Figure 1(c) and the left panel of Figure 2(b) can be directly compared. It can be seen that the proposed cpd_hnf_svd method, even without Newton refinement, can be up to 3 orders of magnitude more accurate (in relative backward error). The differences in accuracy are thus not exclusively attributable to the Newton refinement.

We did not include a visualization of the timings, because they visually look very similar. By adding the Newton refinement, the total execution time over all cases only modestly increased: \(3.4\%\) for eigs and \(3.8\%\) for svd. The relative increase was larger for the small cases while it was close to 0–5% for the large cases (\(n+1 \ge 20\)). Given that the accuracy can be improved by orders of magnitude, we conclude from these experiments that a few steps of Newton refinement of the approximate roots \((\beta _i,\gamma _i)\) is highly recommended.

After deciding that Newton refinement is our default choice and observing in the right panels of Figure 1 there is no significant difference in accuracy between the eigs and svd, we investigate their relative computational performance. In terms of total execution time over all cases, the eigs variant is 29.1% faster. For the largest case this increases to a 40.9% reduction in execution time.

From the above experiments, we conclude that the variant with an iterative eigensolver and Newton refinement is an appropriate default choice. In the remainder of the experiments, cpd_hnf refers to cpd_hnf_eigs+newton, the variant using Arpack’s iterative eigenvalue method applied to the Gram matrix combined with 3 Newton steps to refine the approximate roots.

In the second experiment, we compare cpd_hnf with the current state of the art in direct numerical methods for tensor rank decomposition in terms of accuracy and computational performance. The algorithm developed by Domanov and De Lathauwer [2017];, 2014] is an advanced direct numerical method that shares several high-level characteristics with cpd_hnf:

A Matlab implementation of Domanov and De Lathauwer’s algorithm [Domanov and De Lathauwer 2017] is available as part of the Tensorlab+ repository [Hendrikx et al.] as cpd3_gevd in the domanov2017laa directory. We refer to it as cpd_ddl henceforth. cpd_ddl has an accurate and a fast option for computing the kernel of the auxiliary matrix. We found in our experiments that the dimension of this kernel is often misjudged by the fast implementation. Consequently, it usually increases the size of the auxiliary matrix and causes it to exceed the 16.2 GB memory limit, we put on the size of that matrix. As a concrete statistic, in the configuration with \(d=2\) below, the fast version of cpd_ddl failed in over 75% of the cases. Therefore, we exclusively compare our implementation with the accurate version cpd_ddl.

We generate random rank-\(r\) tensors \(\mathcal {A}\) as described at the start of this section. We apply both cpd_hnf and cpd_ddl to these random tensors for all of the following configurations¹³:and in all cases, we take \(\ell\) \(+\) \(1 = r = \min \lbrace \mathcal {R}(m, n, (d,1)), mn\rbrace\). We do not provide our algorithm with the value of \(d\). For each input tensor, cpd_hnf determines the degree \((d,1)\) or \((1,e)\) automatically as explained in the previous subsection. For each algorithm, we record the relative backward error, total execution time, and the size of the auxiliary matrix whose kernel is computed. As the latter is the dominant operation in both algorithms, its time and memory complexity gives a good estimate of the overall complexity.

For the first configuration cpd_ddl was executed on bigmem nodes of the supercomputer. These are the same as the skylake nodes, except that they are equipped with 768 GB of main memory. The reason was that our Matlab driver routine consumed more than 162 GB of memory while running through all configurations.

The accuracy of cpd_hnf and cpd_ddl in this set of experiments is shown in Figure 2. The newly proposed algorithm is consistently several orders of magnitude more accurate in terms of the relative backward error than the state of the art.

It took over 146 h to obtain the (incomplete) results in the left panel of Figure 2(b). To put this in perspective, we note that the computation for the left panel of Figure 2(a) took a little over 17 h. The missing values inside of the triangles in Figure 2(b) indicate that cpd_ddl wanted to allocate an auxiliary matrix that would require more than 16.2 GB of memory. The same memory constraint was also imposed on cpd_hnf, but here only the largest case with \(d=4\) and \(m=n=14\) could not be treated. The red pixel for \(d = 2, m = 35, n = 13\) in Figure 2(a) indicates that cpd_hnf gave inaccurate results. This is the only case where eigs failed to find a sufficiently accurate nullspace.

The timings of our Julia implementation are shown in Figure 3. As cpd_ddl is implemented in a different programming language, we believe a direct comparison in timings is not opportune. Nevertheless, in both algorithms computing the kernel of the auxiliary matrix has the highest asymptotic time complexity. In cpd_ddl it is an \({\rm L}\) \(\times\) \({\rm L}\) square matrix, and in our algorithm, depending on the choice of degree, \(R_{I}(d,e)\) is an almost-square \({\rm M} \times {\rm N}\) matrix with \({\rm M} \approx {\rm N}\). Therefore, we decided to plot the ratio between the number of elements in the auxiliary matrix of cpd_ddl and cpd_hnf. Figure 4 visualizes this factor \(\mu = {\rm L}^2/({\rm M}{\rm N})\) in a logarithmic scale. This number indicates the fraction of memory that cpd_ddl requires relative to cpd_hnf. Raising it to the power \(\frac{3}{2}\) gives an estimate of the speedup factor in execution time of cpd_hnf relative to cpd_ddl. Based on this estimate and the execution times we logged, it is accurate to state that the newly proposed algorithm outperforms the state of the art by up to two orders of magnitude for larger tensors.

The next experiment illustrates that cpd_hnf can successfully decompose tensors even in the presence of some noise. The setup is as follows. We generate a random rank-\(r\) tensor \(\mathcal {A}\) of size \(150 \times 25 \times 10\) by randomly sampling the \(\alpha _i,\beta _i\), and \(\gamma _i\) as before. Then, we add white Gaussian noise of relative magnitude 10^e for \(e=-1,\ldots ,-15\); that is, we compute \(\mathcal {A}^{\prime } = \mathcal {A}+ 10^{e} \frac{\Vert \mathcal {A}\Vert _F}{\Vert \mathcal {E}\Vert _F} \mathcal {E}\). We provide \(\mathcal {A}^{\prime }\) as input to our algorithm and request a rank-\(r\) decomposition.

The relative backward error between \(\mathcal {A}^{\prime }\) and the computed rank-\(r\) decomposition is shown in a logarithmic scale in Figure 5. Because of our setup, the rank-\(r\) CPD of \(\mathcal {A}\) has relative backward error 10^e. A good tensor decomposition algorithm should thus return a CPD with a backward error of at most 10^e. Remarkably, for tensors with random rank-\(r\) CPDs, the proposed algorithm consistently manages to reach this benchmark when \(e \le -5\). For ranks up to about half the maximum range (from \(r=1\) to 70), it even consistently manages to reach this benchmark for white noise of magnitude at most \(10^{-2}\)! Based on these results, we anticipate that cpd_hnf could be employed as a rank-\(r\) approximation algorithm in the high signal-to-noise regime. We believe this observation warrants further research.

The last experiment illustrates that the reshaping trick combined with a decomposition algorithm that works well in the unbalanced regime (such as cpd_hnf and cpd_ddl [Domanov and De Lathauwer 2014;, 2017]) is a powerful technique for decomposing high-order, high-rank tensors, even with a balanced shape.

As an example, we generated an eighth-order real tensor of size \(7 \times 7 \times 7 \times 7 \times 6 \times 6 \times 5 \times 5\) and rank \(r = 1,000\) with factor matrices whose entries are sampled i.i.d. from a standard normal distribution. Here is what happens:

To our knowledge, this computation represents the first time any tensor decomposition algorithm of any type (i.e., alternating least squares, optimization-based methods, direct algebraic methods, homotopy-based methods, or heuristic methods) successfully decomposes a high-order rank-1000 tensor that cannot be reshaped to an order-3 tensor whose CPD can be computed with a pencil-based algorithm.