Decoupling multivariate polynomials: interconnections between tensorizations

TL;DR

This paper explores the relationships between different tensor-based methods for decoupling multivariate polynomials, revealing their connections and how considering intrinsic tensor structures enhances decomposition uniqueness and applicability.

Contribution

It uncovers the linear relationships between tensor representations in polynomial decoupling and improves decomposition uniqueness by leveraging intrinsic tensor structures.

Findings

Tensors involved are linearly related with interconnected CP decompositions.

Considering intrinsic tensor structure enhances uniqueness of decompositions.

The approach broadens the applicability of tensor-based polynomial decoupling methods.

Abstract

Decoupling multivariate polynomials is useful for obtaining an insight into the workings of a nonlinear mapping, performing parameter reduction, or approximating nonlinear functions. Several different tensor-based approaches have been proposed independently for this task, involving different tensor representations of the functions, and ultimately leading to a canonical polyadic decomposition. We first show that the involved tensors are related by a linear transformation, and that their CP decompositions and uniqueness properties are closely related. This connection provides a way to better assess which of the methods should be favored in certain problem settings, and may be a starting point to unify the two approaches. Second, we show that taking into account the previously ignored intrinsic structure in the tensor decompositions improves the uniqueness properties of the decompositions and thus enlarges the applicability range of the methods.

Equations

$$\mathbf{f}(\mathbf{u})={\left[\begin{array}{ccc}{f}_1(\mathbf{u})& \cdots & f_n(\mathbf{u})\end{array}\right]}^{\top }$$

$$\mathbf{f}(\mathbf{u})=\mathbf{Wg}({\mathbf{V}}^{\top }\mathbf{u}),$$

$$\mathbf{g}(x_1,\dots ,x_r)={\left[\begin{array}{ccc}{g}_1(x_1)& \cdots & g_r(x_r)\end{array}\right]}^{\top },$$

$$g_k(t)=c_{1}t+\cdots +c_d{t}^d.$$

$$\mathbf{f}(\mathbf{u})={\mathbf{w}}_1{g}_1({\mathbf{v}}_1^{\top }\mathbf{u})+\cdots +{\mathbf{w}}_r{g}_r({\mathbf{v}}_r^{\top }\mathbf{u}),$$

$$\begin{array}{rl}{f}_1(u_1,u_2)& =-3u_1^3-9u_1^2{u}_2-27u_1{u}_2^2-15u_2^3-8u_1^2-8u_1{u}_2-20u_2^2+3u_1+9u_2,\\ f_2(u_1,u_2)& =-7u_1^3-6u_1^2{u}_2+6u_1{u}_2^2+7u_2^3+10u_1^2+16u_1{u}_2+10u_2^2-3u_2.\end{array}$$

$$\mathbf{V}=\left[\begin{array}[]{rrr}2&-1&1\\ 1&1&2\end{array}\right],\quad\mbox{and}\quad\mathbf{W}=\left[\begin{array}[]{rrr}0&1&-2\\ -1&0&1\end{array}\right],$$

$$-8u_{1}^{2}-8u_{1}u_{2}-20u_{2}^{2}=\mathbf{u}^{\top}\Psi^{(2)}\mathbf{u},\quad\mbox{where }\Psi^{(2)}=\left[\begin{array}[]{rr}-8&-4\\ -4&-20\end{array}\right].$$

$$p(\mathbf{u})={\Psi }^{(d)}{\bullet }_1\mathbf{u}\cdots {\bullet }_d\mathbf{u}.$$

$$p(u_1,\dots ,u_m)=\sum _{i=1}^r{w}_i(v_{1i}{u}_1+\cdots +v_{mi}{u}_m{)}^d,$$

$${\Psi }^{(d)}=\sum _{i=1}^r{w}_i({\mathbf{v}}_i\circ \cdots \circ {\mathbf{v}}_i).$$

$$\left[\begin{array}[]{rr}-8&-4\\ -4&-20\end{array}\right]=\left[\begin{array}[]{rr}-2&2\\ 2&4\end{array}\right]\left[\begin{array}[]{rr}-1&0\\ 0&-1\end{array}\right]\left[\begin{array}[]{rr}-2&2\\ 2&4\end{array}\right],$$

$$p(u_1,u_2)=-(-2u_1+2u_2{)}^2-(2u_1+4u_2{)}^2.$$

$$\left[\begin{array}[]{rr}-8&-4\\ -4&-20\end{array}\right]\approx\left[\begin{array}[]{rr}0.2898&-0.9571\\ 0.9571&0.2898\end{array}\right]\left[\begin{array}[]{rr}-21.2111&0\\ 0&-6.7889\end{array}\right]\left[\begin{array}[]{rr}0.2898&0.9571\\ -0.9571&0.2898\end{array}\right]$$

$$\begin{array}[]{rcl}p_{1}(u_{1},\ldots,u_{m})&=&\Psi^{(d)}_{1}\bullet_{1}\mathbf{u}\cdots\bullet_{d}\mathbf{u},\\ &\vdots&\\ p_{n}(u_{1},\ldots,u_{m})&=&\Psi^{(d)}_{n}\bullet_{1}\mathbf{u}\cdots\bullet_{d}\mathbf{u}.\\ \end{array}$$

$$p(\mathbf{u})={\mathbf{u}}^{\top }{\Psi }^{(1)}+{\mathbf{u}}^{\top }{\Psi }^{(2)}\mathbf{u}+{\Psi }^{(3)}{\bullet }_1\mathbf{u}{\bullet }_2\mathbf{u}{\bullet }_3\mathbf{u}+\cdots +{\Psi }^{(d)}{\bullet }_1\mathbf{u}\cdots {\bullet }_d\mathbf{u},$$

$$f_{1}(u_{1},u_{2})=\mathbf{u}^{\top}\left[\begin{array}[]{r}3\\ 9\end{array}\right]+\mathbf{u}^{\top}\left[\begin{array}[]{rr}-8&-4\\ -4&-20\end{array}\right]\mathbf{u}+\Psi^{(3)}\bullet_{1}\mathbf{u}\bullet_{2}\mathbf{u}\bullet_{3}\mathbf{u},$$

$$\Psi^{(3)}_{:,:,1}=\left[\begin{array}[]{rr}-3&-3\\ -3&-9\end{array}\right],\quad\quad\quad\Psi^{(3)}_{:,:,2}=\left[\begin{array}[]{rr}-3&-9\\ -9&-15\end{array}\right],$$

$$f_{2}(u_{1},u_{2})=\mathbf{u}^{\top}\left[\begin{array}[]{r}0\\ -3\end{array}\right]+\mathbf{u}^{\top}\left[\begin{array}[]{rr}10&8\\ 8&10\end{array}\right]\mathbf{u}+\Psi^{(3)}\bullet_{1}\mathbf{u}\bullet_{2}\mathbf{u}\bullet_{3}\mathbf{u},$$

$$\Psi^{(3)}_{:,:,1}=\left[\begin{array}[]{rr}-7&-2\\ -2&2\end{array}\right],\quad\quad\quad\Psi^{(3)}_{:,:,2}=\left[\begin{array}[]{rr}-2&2\\ 2&7\end{array}\right].$$

$$\begin{array}[]{rcl}p_{1}(u_{1},\ldots,u_{m})&=&\mathbf{u}^{\top}\Psi^{(1)}_{1}+\cdots+\Psi^{(d)}_{1}\bullet_{1}\mathbf{u}\cdots\bullet_{d}\mathbf{u},\\ &\vdots&\\ p_{n}(u_{1},\ldots,u_{m})&=&\mathbf{u}^{\top}\Psi^{(1)}_{n}+\cdots+\Psi^{(d)}_{n}\bullet_{1}\mathbf{u}\cdots\bullet_{d}\mathbf{u},\end{array}$$

$${\mathcal{Q}}_{i,:,:}:=\Psi (f_i),$$

$$\Psi(p)=\left[\begin{array}[]{c|c|c|c|c}\Psi^{(1)}&\Psi^{(2)}&\Psi^{(3)}_{(1)}&\cdots&\Psi^{(d)}_{(1)}\end{array}\right],$$

$$p(u_1,u_2)=a_1{u}_1+a_2{u}_2+b_1{u}_1^2+2b_2{u}_1{u}_2+b3u_2^2+d_1{u}_1^3+3d_2{u}_1^2{u}_2+3d_3{u}_1{u}_2^2+d_4{u}_3^3$$

$$p(u_1,u_2)={\mathbf{u}}^{\top }\left[\begin{array}{c}{a}_1\\ a_2\end{array}\right]+{\mathbf{u}}^{\top }\left[\begin{array}{cc}{b}_1& b_2\\ b_2& b_3\end{array}\right]\mathbf{u}+{\Psi }^{(3)}{\bullet }_1\mathbf{u}{\bullet }_2\mathbf{u}{\bullet }_3\mathbf{u},$$

$${\Psi }_{:,:,1}^{(3)}=\left[\begin{array}{cc}{d}_1& d_2\\ d_2& d_3\end{array}\right],{\Psi }_{:,:,2}^{(3)}=\left[\begin{array}{cc}{d}_2& d_3\\ d_3& d_4\end{array}\right].$$

$$\Psi(p)=\left[\begin{array}[]{c|cc|cccc}a_{1}&b_{1}&b_{2}&d_{1}&d_{2}&d_{2}&d_{3}\\ a_{2}&b_{2}&b_{3}&d_{2}&d_{3}&d_{3}&d_{4}\end{array}\right].$$

$$\mathcal{Q}_{1,:,:}=\Psi(f_{1})=\left[\begin{array}[]{c|cc|cccc}3&-8&-4&-3&-3&-3&-9\\ 9&-4&-20&-3&-9&-9&-15\end{array}\right],$$

$$\mathcal{Q}_{2,:,:}=\Psi(f_{2})=\left[\begin{array}[]{c|cc|cccc}0&10&8&-7&-2&-2&2\\ -3&8&10&-2&2&2&7\end{array}\right].$$

$$\mathcal{Q}=\sum _{k=1}^r{\mathbf{w}}_k\circ {\mathbf{v}}_k\circ {\mathbf{z}}_k,$$

Full text

We first show that the involved tensors are related by a linear transformation, and that their CP decompositions and uniqueness properties are closely related. This connection provides a way to better assess which of the methods should be favored in certain problem settings, and may be a starting point to unify the two approaches. Second, we show that taking into account the previously ignored intrinsic structure in the tensor decompositions improves the uniqueness properties of the decompositions and thus enlarges the applicability range of the methods.

1 Introduction

Representing a nonlinear function in a simpler way can provide an insight into its inner workings, reduce the parametric complexity, or facilitate function approximation. One of the successful examples are tensor decompositions, such as the canonical polyadic (CP) decomposition that can be viewed as a decomposition of functions into a sum of separable functions [1]. Tensor decompositions found many applications in signal/image processing, chemometrics, physics, machine learning, to name a few [2, 3, 4, 5, 6]. In these applications, tensors either appear naturally due to multi-dimensionality of data [7], or the data can be tensorized, i.e., a higher-order tensor is constructed from data [8].

In this paper, we focus on the task of decoupling a set of polynomial vector functions, that is, decomposing a set of multivariate real polynomials into linear combinations of univariate polynomials in linear forms of the input variables. This task has attracted a spark of research attention over the last years, motivated by several applications, such as system identification [9, 10, 11, 12, 13, 14], approximation theory [15, 16, 17], and neural networks [18]. Restricting polynomial decoupling to a single homogeneous polynomial is equivalent to the well-known Waring decomposition [19, 20], but some generalizations to non-homogeneous polynomials or joint Waring decompositions are studied as well [21, 22] and [23, 24, 11].

Several tensor-based approaches were proposed for computing a decoupled representation of a given function [12, 13, 25, 26, 24]. These solutions can be categorized into two classes. The methods [12, 13, 25, 26] build a tensor from the polynomial coefficients, whereas the method of [24] builds a tensor from the Jacobian matrices of the functions, evaluated at a set of sampling points. Ultimately, all methods boil down to a canonical polyadic decomposition (CP decomposition) of the constructed tensor to retrieve a decoupled representation in which the nonlinearities occur as univariate polynomial mappings.

The benefit of using a tensor-based approach for decoupling is twofold. First, ‘tensorization’ procedures often lead to (essentially) uniquely decomposable tensors [8], i.e., ensuring that identifiable structures can be retrieved. Second, by solving the decoupling problem as a CP decomposition, one can use recent widely available and robust numerical tools, such as Tensorlab for MATLAB [27] (or alternatives [28, 29]).

This paper specifically focuses on the two tensorization methods [26] and [24]. Although both associated tensors have a particular structure, both approaches seem quite different in nature, and each of the methods has distinct advantages over the other one. For instance, the coefficient-based methods [12, 13, 25, 26] require several high-order tensors (or their matricizations) for polynomials of high degrees, whereas [24] involves a single third-order tensor only. Coefficient-based approaches can easily deal with single polynomials, whereas [24] would in that case not be able to take advantage of the uniqueness properties of the CP decomposition, as the tensor of Jacobian matrices is then a matrix composed of gradient vectors. On the other hand, the approach of [24] can be applied to non-polynomial functions, which may in some cases be of interest, e.g., in [9] a neural network was decoupled.

We aim at obtaining a deeper understanding of the connections between the solution approaches. This is profitable when extending the applicability range of the methods, e.g., when moving from polynomials to any differentiable functions. Furthermore, such connections may provide a way to transfer theoretical properties from one formulation to another. For example, as we argue in Section 6, exploring the previously ignored structure in the tensor decomposition in one of the settings enlarges the range of decomposable functions. This knowledge may lead to improved algorithms in another setting as well.

The remainder of this article is organized as follows: Section 2 formalizes the problem of decoupling multivariate polynomials. Section 3 explains the link between the decoupling problem and the symmetric tensor decomposition problem. Section 4 discusses the construction of the tensor of unfoldings [26] and the Jacobian tensor [24]. Section 5 presents our first contribution, namely the relation between the two tensorizations. The second main contribution of the paper is Section 6, which clarifies the need of dealing with structure in the decompositions and proposes a coupled CP decomposition approach for solving the structured problem. Section 7 draws the conclusions and points out open problems for future work.

Notation

Scalars are denoted by lowercase or uppercase letters. Vectors are denoted by lowercase boldface letters, e.g., $\mathbf{u}$. Elements of a vector are denoted by lowercase letters with an index as subscript, e.g., \mathbf{x}=\left[\begin{array}[]{ccc}x_{1}&\cdots&x_{m}\end{array}\right]^{\top}. Matrices are denoted by uppercase boldface letters, e.g., $\mathbf{V}$. The entry in the $i$-th row and $j$-th column of a matrix $\mathbf{V}$ is denoted by $v_{ij}$, and the matrix $\mathbf{V}\in\mathbb{R}^{m\times r}$ may be represented by its columns \mathbf{V}=\left[\begin{array}[]{ccc}\mathbf{v}_{1}&\cdots&\mathbf{v}_{r}\end{array}\right]. The Kronecker product of matrices is denoted by “$\otimes$”.

Tensors of order $d$ are denoted by uppercase caligraphical letters, e.g., $\mathcal{J}\in\mathbb{R}^{n\times m\times N}$. The outer product is denoted by “$\circ$” and is defined as follows: For $\mathcal{T}=\mathbf{u}\circ\mathbf{v}\circ\mathbf{w}$, the entry in position $(i,j,k)$ is equal to $u_{i}v_{j}w_{k}$. The canonical polyadic (CP) decomposition expresses a tensor $\mathcal{T}$ as a (minimal) sum of rank-one tensor terms [30, 31, 2] as $\mathcal{T}=\sum_{i=1}^{R}\mathbf{u}_{i}\circ\mathbf{v}_{i}\circ\mathbf{w}_{i}$, and is sometimes denoted in a short-hand notation as $\mathcal{T}=\llbracket\mathbf{U},\mathbf{V},\mathbf{W}\rrbracket$. The CP rank $r$ is defined as the (minimal) number of terms that is required to represent $\mathcal{T}$ as a sum of $r$ rank-one terms. To refer to elements of matrices or tensors, or subsets thereof, we may use MATLAB-like index notation (including MATLAB’s colon wildcard): for instance, $\mathcal{T}_{i,j,k,\ell}$ is the element at position $(i,j,k,\ell)$ of a fourth-order tensor $\mathcal{T}$, and $\mathcal{T}_{:,:,2}$ is the second frontal slice of a third-order tensor $\mathcal{T}$. The mode-$n$ product is denoted by “$\bullet_{n}$” and is defined as follows. Let $\mathcal{X}$ be an $I_{1}\times I_{2}\times\cdots\times I_{N}$ tensor, and let $\boldsymbol{u}$ be a vector of length $I_{n}$, then we have $\left(\mathcal{X}\bullet_{n}\boldsymbol{u}^{\top}\right)_{i_{1}\cdots i_{n-1}i_{n+1}\cdots i_{N}}=\sum_{i_{n}=1}^{I_{n}}x_{i_{1}i_{2}\cdots i_{N}}u_{i_{n}}$. Notice that the result is a tensor of order $N-1$, as mode $n$ is summed out. Similarly, for an $I_{1}\times I_{2}\times\cdots\times I_{N}$ tensor $\mathcal{X}$ and a matrix $\mathbf{M}\in J\times I_{n}$, the mode-$n$ product is defined as $\left(\mathcal{X}\bullet_{n}\mathbf{M}^{\top}\right)_{i_{1}\cdots i_{n-1}ji_{n+1}\cdots i_{N}}=\sum_{i_{n}=1}^{I_{n}}x_{i_{1}i_{2}\cdots i_{N}}m_{j,i_{n}}$. Let $\mathop{\operator@font vec}(\mathcal{T})$ denote the column-major vectorization of a tensor $\mathcal{T}$. The first-mode unfolding of an $I_{1}\times I_{2}\times\cdots\times I_{N}$ tensor $\mathcal{X}$ is the matrix $\mathcal{X}_{(1)}$ of size $I_{1}\times I_{2}\cdots I_{N}$, where each row is the vectorized slice of the tensor $\mathcal{X}$, i.e. $(\mathcal{X}_{(1)})_{i,:}=\left(\mathop{\operator@font vec}(\mathcal{X}_{i,:,\ldots,:})\right)^{\top}$ (see, for example, [2] for more details).

2 The polynomial decoupling model

First, we describe the model, following the notation of [24] as illustrated in Fig. 1. Consider a multivariate polynomial map $\mathbf{f}:\mathbb{R}^{m}\to\mathbb{R}^{n}$, i.e., a vector

[TABLE]

of multivariate polynomials (of total degree at most $d$) in variables $\mathbf{u}=\begin{bmatrix}u_{1}&\cdots&u_{m}\end{bmatrix}^{\top}$. We say that $\mathbf{f}$ has a decoupled representation, if it can be expressed as

[TABLE]

where $\mathbf{V}\in\mathbb{R}^{m\times r},\mathbf{W}\in\mathbb{R}^{n\times r}$ are transformation matrices, and $\mathbf{g}:\mathbb{R}^{r}\to\mathbb{R}^{r}$ is defined as

[TABLE]

where $g_{k}:\mathbb{R}\to\mathbb{R}$ are univariate polynomials of degree at most $d$, i.e.,

[TABLE]

Note that we omitted the constant terms of the polynomials, since they are not uniquely identifiable [24]. In this paper we limit ourselves to the model (1).

The decoupled representation (1) can be also equivalently rewritten as

[TABLE]

where $\mathbf{v}_{k}$ and $\mathbf{w}_{k}$ are the columns of $\mathbf{V}$ and $\mathbf{W}$, respectively. As shown in [32, 33], the decomposition (3) is a special case of the $X$-rank decomposition [34, §5.2.1], where the set of “rank-one” terms is the set of polynomial maps of the form $\mathbf{w}g(\mathbf{v}^{\top}\mathbf{u})$. The $X$-rank framework is useful [33] for studying the identifiability of the model (3).

The following example shows a decoupled representation for a simple case. This example will be used throughout the paper to illustrate the main ideas of the various aspects that we will explore.

Example 1**.**

Consider a function \mathbf{f}(\mathbf{u})=\left[\begin{array}[]{cc}f_{1}(u_{1},u_{2})&f_{2}(u_{1},u_{2})\end{array}\right]^{\top} given as

[TABLE]

It can be verified that $\mathbf{f}$ has a decomposition (3) with $m=n=2$ and $r=3$ as

[TABLE]

*and $g_{1}(x_{1})=x_{1}^{3}-2x_{1}^{2}-x_{1}$, $g_{2}(x_{2})=x_{2}^{3}-4x_{2}^{2}+x_{2}$, $g_{3}(x_{3})=x_{3}^{3}+2x_{3}^{2}-2x_{3}$ (see Figure 2). *

Remark 1**.**

In general, the coupled representation $\mathbf{f}(\mathbf{u})$ has $n\left({m+d\choose d}-1\right)$ coefficients, while the decoupled representation $\mathbf{W}\mathbf{g}(\mathbf{V}^{\top}\mathbf{u})$ has $r(m+n+d)$ coefficients. Due to the combinatorial increase of the number of coefficients in the coupled representation, the decoupled representation is especially beneficial for large values of $m$, $n$, and $d$. But even for small values of $m$, $n$, and $d$, the parametric reduction can be significant, for example, if $m=n=3$, $d=5$, and $r=3$, the coupled representation has 168 coefficients, while the decoupled one has only 36 coefficients.

3 Decoupling polynomials and symmetric tensor decompositions

Let us review some well-known facts that connect polynomials with symmetric tensors [35, 36], and that connect some special cases of the representation (1) with symmetric tensor decompositions.

3.1 Homogeneous polynomials, symmetric tensors and Waring decomposition

It is well-known that there is a one-to-one correspondence between homogeneous polynomials and symmetric tensors [35]. For instance, the polynomial $-8u_{1}^{2}-8u_{1}u_{2}-20u_{2}^{2}$ can be written as

[TABLE]

In general, let $p^{(d)}(u_{1},\ldots,u_{m})$ be a homogeneous polynomial (also called a $d$-ary form) of degree $d$ in $m$ variables. Then there is a unique symmetric tensor $\Psi^{(d)}$ of order $d$ and dimension $m$ such that

[TABLE]

Next, it is easy to see that the decoupling problem for the polynomial (5) takes the form

[TABLE]

which is known as the Waring decomposition [19, 20] of $p(u_{1},\ldots,u_{m})$. The Waring decomposition, in its turn, is equivalent to the symmetric CP decomposition of $\Psi^{(d)}$:

[TABLE]

The symmetric CP decomposition of $\Psi^{(d)}$ reveals possible values for the unknowns $v_{ij}$ and $w_{i}$.

Example 2**.**

Consider the polynomial given in (4). Then the corresponding symmetric matrix $\Psi^{(2)}$ admits the decomposition

[TABLE]

such that $p(u_{1},u_{2})=\mathbf{u}^{\top}\Psi^{(2)}\mathbf{u}$ has the Waring decomposition

[TABLE]

Notice that the symmetric decomposition of $\Psi^{(2)}$ from Example 2 is not unique (nor ‘essentially unique’ [2]). Indeed, the eigenvalue decomposition

[TABLE]

provides another valid factorization. For $d>2$, however, the Waring decomposition (6) possesses uniqueness properties even in the case of quite large ranks [37, 38].

Along the same lines, it is possible to decouple jointly several homogeneous polynomials. Consider the case of $n$ homogeneous polynomials of degree $d$, denoted by

[TABLE]

Then the decoupling problem (1) corresponds to the simultaneous Waring decomposition of several forms or, equivalently, the coupled CP decomposition of several symmetric tensors. The rank and identifiability properties of simultaneous Waring decompositions were also studied in the literature, see [23, 37, 39] and references therein.

3.2 The case of non-homogeneous polynomials

Next, consider the case of non-homogeneous polynomials. Any non-homogeneous polynomial of degree $d$ can hence be written as

[TABLE]

where $\Psi^{(1)}\in\mathbb{R}^{m}$ , $\Psi^{(2)}\in\mathbb{R}^{m\times m}$ is a symmetric matrix, and each $\Psi^{(s)}\in\mathbb{R}^{m\times\cdots\times m}$, $3\leq s\leq d$, is a symmetric tensor of order $s$.

Example 3**.**

We continue Example 1. We can write $f_{1}(u_{1},u_{2})$ and $f_{2}(u_{1},u_{2})$ as

[TABLE]

with

[TABLE]

and

[TABLE]

with

[TABLE]

The decomposition of a single non-homogeneous polynomial as in (3) is hence equivalent to joint decomposition of several symmetric tensors but of different orders [32].

Finally, several non-homogeneous polynomials can be jointly decomposed in a similar way. Consider $n$ non-homogeneous polynomials of maximal degree $d$, denoted as

[TABLE]

The full decomposition in (1) can be also viewed as a coupled tensor decomposition, which will be presented in Section 6.2.

4 Tensorizations and their decompositions

In this section, we recall tensorizations proposed in the literature to find the decomposition (1) by a CP decomposition of a single tensor constructed from $\mathbf{f}$, namely the tensorizations of [26] and [24]. We recall basic properties and give short proofs for completeness, although these proofs are already present in [26, 24]. We also use a slightly different notation to simplify the exposition.

4.1 Tensor of unfoldings [26]

The above link between polynomials, (partially) symmetric tensors and their CP decompositions gives rise to the tensorization approach of [26], in which a tensor is constructed from the coefficients of the polynomials $f_{1}(u_{1},\ldots,u_{m})$ up to $f_{n}(u_{1},\ldots,u_{m})$. This tensorization offers the advantage that several polynomials can be represented as a single tensor, and the decoupling task can be solved using a single (but structured) CP decomposition. In this approach, the tensor (shown in Figure 3) is constructed from the coefficients of the polynomial map of degree $d$, as follows:

The tensor has size $n\times m\times\delta$, where $\delta=\sum\limits_{k=1}^{d}m^{k-1}$.

• 2.

The tensor is constructed by slices

[TABLE]

where $\Psi$ is a structured $m\times\delta$ matrix built from the coefficients of $f_{i}(\mathbf{u})$.

Now let us describe the construction of the structured coefficient matrix $\Psi(p)$ for a given polynomial of degree $d$. Recall that each such polynomial can be written as in (9), where $\Psi^{(1)}\in\mathbb{R}^{m}$ , $\Psi^{(2)}\in\mathbb{R}^{m\times m}$ is a symmetric matrix and $\Psi^{(s)}\in\mathbb{R}^{m\times\cdots\times m}$ are symmetric tensors of order $s$. Then the matrix $\Psi(p)\in\mathbb{R}^{m\times\delta}$ is constructed111In [26] the linear term is skipped, and $\delta=\sum\limits_{k=2}^{d}m^{d-1}$. In [40] the matrix $\Psi$ is denoted as $\Gamma$. as

[TABLE]

where $\mathcal{G}_{(1)}$ denotes the first-mode unfolding of a tensor $\mathcal{G}$.

Example 4**.**

A third-degree polynomial in two variables

[TABLE]

has the representation

[TABLE]

where

[TABLE]

By putting all the unfoldings together, we get

[TABLE]

Hence, for $f_{1}$ and $f_{2}$ in Example 1, the slices of the tensor $\mathcal{Q}$ are given by

[TABLE]

and

[TABLE]

As proved in [26], the tensor $\mathcal{Q}$ has a CP decomposition, which reveals the decomposition (1). We repeat here a simplified version of the proof for completeness.

Lemma 1**.**

For the polynomial map (1), the tensor $\mathcal{Q}$ has the following CP decomposition:

[TABLE]

where

[TABLE]

Proof.

Consider $q_{k}(\mathbf{u}):=g_{k}(\mathbf{v}_{k}^{\top}\mathbf{u})$, where $g_{k}$ is as in (2). Easy calculations show that

[TABLE]

see also [40, eqn. (A.7)]. Since, from (3) $f_{i}(\mathbf{u})=\sum\limits_{k=1}^{r}(\mathbf{w}_{k})_{i}q_{k}(\mathbf{u})$, we have that

[TABLE]

which implies (14). ∎

Example 5**.**

We continue Examples 1, 4. The Kronecker products of columns of $\mathbf{V}$ are:

[TABLE]

Hence, the matrix $\mathbf{Z}=\begin{bmatrix}\mathbf{z}_{1}&\mathbf{z}_{2}&\mathbf{z}_{3}\end{bmatrix}$ is given by

[TABLE]

4.2 The tensor of Jacobian matrices of [24]

The tensorization method of [24] does not use the coefficients of $\mathbf{f}(\mathbf{u})$ directly, but proceeds by collecting the first-order information of $\mathbf{f}(\mathbf{u})$ (i.e., the partial derivatives) in a set of sampling points. The thusly obtained Jacobian matrices are arranged into a third-order tensor, of which the CP decomposition reveals the decomposition (1).

As in [24], we consider the Jacobian of $\mathbf{f}$:

[TABLE]

Using Lemma 2, the tensorization is constructed as follows (see Figure 4):

$N$ points $\mathbf{u}^{(1)},\ldots,\mathbf{u}^{(N)}\in\mathbb{R}^{m}$ are chosen (so-called sampling points).

• 2.

An $n\times m\times N$ tensor $\mathcal{J}$ is constructed by stacking the Jacobian evaluations at $\mathbf{u}^{(k)}$

[TABLE]

Example 6**.**

*We continue Example 1. As a set of sampling points, we choose *

[TABLE]

By evaluating $\mathbf{J}_{\mathbf{f}}(\mathbf{u})$ at these points, we get the tensor $\mathcal{J}$ given by

[TABLE]

If $\mathbf{f}(\mathbf{u})$ has a decoupled representation (1), the following lemma holds true.

Lemma 2** ([24, Lemma 2.1]).**

The first order derivatives of (1) are given by

[TABLE]

where $g_{i}^{\prime}(t):=\frac{dg_{i}}{dt}(t)$.

The proof, given in [24], follows by chain rule:

[TABLE]

By Lemma 2, the evaluations of the Jacobians can be jointly factorized:

[TABLE]

where $\mathbf{D}^{(k)}=\mathop{\operator@font diag}(g^{\prime}_{1}(\mathbf{v}_{1}^{\top}\mathbf{u}^{(k)}),\ldots,g^{\prime}_{r}(\mathbf{v}_{r}^{\top}\mathbf{u}^{(k)}))$. Therefore, $\mathcal{J}$ admits a CP decomposition

[TABLE]

where $\mathbf{w}_{k}$, $\mathbf{v}_{k}$ are as in (3), and $\mathbf{h}_{k}$ contains the evaluations of $g^{\prime}_{k}(\mathbf{v}_{k}^{\top}\mathbf{u})$ in $\mathbf{u}^{(1)},\ldots,\mathbf{u}^{(N)}$:

[TABLE]

Example 7**.**

We continue Examples 1 and 6. By differentiation, we get

[TABLE]

and hence, by substitution,

[TABLE]

Straightforward calculations show indeed that $\mathcal{J}$ given in (17) admits a decomposition (20) with $\mathbf{H}$ as in (22).

5 Relation between tensorizations $\mathcal{J}$ and $\mathcal{Q}$

In this section, we show how CP decompositions of (14) and (20) are related. Moreover, we establish the relation between the ranks of the tensors and uniqueness of CP decompositions.

First, we show the relation between the vectors $\mathbf{z}_{k}$ and $\mathbf{h}_{k}$, defined in (15) and (21), respectively.

We give the proof of this basic fact for completeness.

Lemma 3**.**

*The vectors $\mathbf{z}_{k}$ and $\mathbf{h}_{k}$ defined in (15) and (21), respectively, satisfy *

[TABLE]

where $\mathbf{A}\in\mathbb{R}^{\delta\times N}$ is a Vandermonde-like matrix whose columns are

[TABLE]

Proof.

Recall that by the properties of the Kronecker product

[TABLE]

Then from (15) have that

[TABLE]

where the last equality follows from (21) and the fact that

[TABLE]

∎

Example 8**.**

In Example 6, the matrix $\mathbf{A}$ can be found as

[TABLE]

It is easy to see that $\mathbf{H}=\mathbf{A}^{\top}\mathbf{Z}$.

As a consequence, we get that the two tensors and their ranks are also related.

Theorem 1**.**

For any polynomial map $\mathbf{f}$, $\mathcal{J}$ and $\mathcal{Q}$ are related as

[TABLE] 2. 2.

The rank of $\mathbf{A}$ is bounded as

[TABLE]

In addition, if $M\leq N$, and $M$ points in $\{\mathbf{u}^{(j)}\}$ are in general position, then $\mathop{\operator@font rank}\mathbf{A}=M$. For example, if points $\{\mathbf{u}^{(j)}\}$ are independent and sampled from a continuous probability distribution, then $\mathop{\operator@font rank}\mathbf{A}=M$ with probability $1$. 3. 3.

If $\mathbf{A}$ has maximal possible rank (i.e. $\mathop{\operator@font rank}\mathbf{A}=M$), then

[TABLE]

and all the minimal CP decompositions differ only by the third factors, which are linked as in (23). Moreover, if the CP decomposition of $\mathcal{Q}$ is unique, then the CP decomposition of $\mathcal{J}$ is also unique.

Proof of Lemma 3.

Let us express $g^{\prime}_{k}(\mathbf{v}_{k}^{\top}\mathbf{u})$ in an explicit form. First, $g^{\prime}_{k}(t)=c_{k,1}+2c_{k,2}t+3c_{k,3}t^{2}+\cdots+dc_{k,d}t^{d-1}$, from which it follows that

[TABLE]

Since $(\mathbf{v}^{\top}\mathbf{u})^{s}=(\mathbf{v}\otimes\cdots\otimes\mathbf{v})^{\top}(\mathbf{u}\otimes\cdots\otimes\mathbf{u})$, the $j$-th element of $\mathbf{h}_{k}$ is equal to

[TABLE]

which completes the proof. ∎

Proof of Theorem 1.

1. First, any polynomial map $\mathbf{f}$ can be decomposed as (1) with $r$ sufficiently large. Let us take such a decomposition; then it holds that

[TABLE]

where the last equality follows from (23).

2. By construction, each element in the image of $\mathbf{A}$ lies in the following subspace:

[TABLE]

Taking into account that the dimension of the space of $m\times\cdots\times m$ symmetric tensors of order $s$ is $\binom{m+s-1}{s}$, we get that the maximal possible rank of $\mathbf{A}$ is

[TABLE]

Next, from (24), we have that the $k$-th column contains evaluations of all $M$ monomials $\{u^{j_{1}}_{1}\cdots u^{j_{m}}_{m}\leq d-1\}^{j_{1}+\cdots+j_{m}\leq d}_{j_{1},\ldots,j_{m}=0}$ at a point $\mathbf{u}^{(k)}$ (scaled by a constant). If, without loss of generality, the first $M$ points $\{\mathbf{u}^{(k)}\}_{k=1}^{M}$ are in general position, then the columns of $\mathbf{A}$ corresponding to different monomials are linearly independent by [41, Multiplicity One Theorem], hence $\mathop{\operator@font rank}\mathbf{A}=M$.

Note that each tube $\mathcal{Q}_{i,j,:}$ of the tensor $\mathcal{Q}$, by construction, lies in $\mathscr{A}$. If $\mathop{\operator@font rank}\mathbf{A}=M$, then its row span coincides with $\mathscr{A}$. Hence the following identity holds true:

[TABLE]

The remaining properties follow from (27) and (25). ∎

6 Structured tensor decompositions

6.1 From CPD to a decomposition with structured rank-one terms

The CP decomposition of $\mathcal{J}$ and $\mathcal{Q}$, although related, are not always equivalent to the original decomposition (3). This happens because there are still nontrivial linear dependencies between the elements of $\mathcal{Q}$ and $\mathcal{J}$. In what follows, we establish relationships between the CP decompositions and the original decomposition (3).

First, we prove that for the rank-one case, these decompositions coincide.

Proposition 1**.**

Consider a polynomial map $\mathbf{f}(\mathbf{u})$ of degree $d$, and the tensor $\mathcal{Q}$ built from it. Then the following holds

[TABLE]

where $\mathbf{w}\in\mathbb{R}^{n}$, $\mathbf{v}\in\mathbb{R}^{m}$ and $g(t)$ is a polynomial of degree $d$.

Proof.

The $\boxed{\Leftarrow}$ follows from Lemma 1. Let us prove the $\boxed{\Rightarrow}$ part. Assume that

[TABLE]

First, since the tensor $\mathcal{Q}$ contains all the coefficients of the derivatives, we have that there exists a polynomial $\widetilde{f}(\mathbf{u})$ such that $\nabla f_{k}(\mathbf{u})=(\mathbf{w})_{k}\nabla\widetilde{f}(\mathbf{u})$. Since the polynomials $f_{k}(\mathbf{u})$ do not have constant terms, we have that

[TABLE]

where $\Psi(\widetilde{f})=\mathbf{v}\mathbf{y}^{\top}$.

Next, let us show that the polynomial $\widetilde{f}$ should necessarily the form $\widetilde{f}(\mathbf{u})={g}(\mathbf{v}^{\top}\mathbf{u})$. Since $\Psi(\widetilde{f})=\mathbf{v}\mathbf{y}^{\top}$, then it follows from (11) that all the unfoldings $\Psi^{(1)}$, $\Psi^{(2)}$, $\Psi^{(3)}_{(1)},\cdots,\Psi^{(d)}_{(1)}$ have rank at most one and their column space is spanned by the vector $\mathbf{v}$. Therefore, we have that

[TABLE]

and hence $\widetilde{f}(\mathbf{u})={g}(\mathbf{v}^{\top}\mathbf{u})$ where

[TABLE]

which completes the proof. ∎

Remark 2**.**

The fact that $\mathop{\operator@font rank}\Psi(\widetilde{f})\leq 1$ implies $\widetilde{f}(\mathbf{u})=g(\mathbf{v}^{\top}\mathbf{u})$ also can be proved alternatively, by noting that the matrix $\Psi{\widetilde{f}}$, after removing duplicate columns, can be reduced to the form $S(f)$ in [33, Proposition 22]. Hence, by [33, Proposition 4.1], the polynomial $\widetilde{f}$ has necessarily the form $\widetilde{f}(\mathbf{u})=g(\mathbf{v}^{\top}\mathbf{u})$. However, this alternative proof requires introducing extra notation, which would be much longer that the proof presented in this paper.

Corollary 1**.**

If the $N$ sampling points are chosen such that the $\mathop{\operator@font rank}(\mathbf{A})=M$, then

[TABLE]

As a corollary of Proposition 1, we get that the original polynomial decomposition (3) is equivalent to a structured CP decomposition.

Corollary 2**.**

Let $\mathcal{L}_{\mathcal{Q}}\subset\mathbb{R}^{n\times m\times\delta}$ be the linear subspace of tensors with the structure of $\mathcal{Q}$. Let the sampling points be chosen such that $\mathop{\operator@font rank}\mathbf{A}=M$, and $\mathcal{L}_{\mathcal{J}}\subset\mathbb{R}^{n\times m\times N}$ be the linear subspace of tensors with the structure of $\mathcal{J}$.

Then the following three statements are equivalent:

the polynomial map $\mathbf{f}(\mathbf{u})$ admits a decomposition (3); 2. 2.

the tensor $\mathcal{Q}(\mathbf{f})$ admits the structured CP decomposition

[TABLE] 3. 3.

the tensor $\mathcal{J}(\mathbf{f})$ admits the structured CP decomposition

[TABLE]

The structure constraint is important: indeed, the CP decomposition of the tensor $\mathcal{Q}$ or $\mathcal{J}$ is not necessarily structured. In general, we do not know even if the CP rank is equal to the structured CP rank (minimal number of terms in (28) or (29)). This is similar to the Comon’s conjecture [35, §5] about symmetric tensors: it is not known whether the symmetric rank of a symmetric tensor equals its non-symmetric rank.

However, if the CP decomposition of a tensor is unique (for example, if it satisfies Kruskal’s uniqueness conditions), then it should necessarily be a structured CP decomposition.

6.2 Computing coupled/structured CP decomposition

Earlier attempts to tackle the structured case were made by [25, 13] and [42, §8, pp. 133–136]. The attempts of [25, 13] have the disadvantage that a tensor is built that has missing values, which increase in number as the polynomial degree grows. The attempt of [42] consisted of parameterizing the internal nonlinear functions $g_{k}$ using their coefficients. Although this seems a promising approach, it turned out to be problematic in practice to build a working algorithm, as the decoupling method led to strongly nonlinear/nonconvex optimization problems.

We propose to tackle the problem by solving a coupled and structured CP decomposition instead. First, let us consider simultaneous decoupling of homogeneous polynomials (8). Let us arrange the $\Psi^{(d)}_{i}$, for $i=1,\ldots,n$ into a tensor $\mathcal{T}^{d}$, such that $\mathcal{T}^{d}_{i,:,\ldots,:}=\Psi^{(d)}_{i}$, for all $i=1,\ldots,n$. Then it is easy to verify that $\mathcal{T}^{d}$ admits a partially symmetric CP decomposition

[TABLE]

which, in our decoupled representation (1), takes the form $\mathbf{W}\mathbf{g}(\mathbf{V}^{\top}\mathbf{u})$, where \mathbf{g}(\mathbf{x})=\left[\begin{array}[]{ccc}x_{1}^{d}&\cdots&x_{r}^{d}\end{array}\right]^{\top}.

Decoupling non-homogeneous polynomials can be achieved by means of a coupled structured CP decomposition of the $\mathcal{T}^{d}$ tensors. Let us arrange all $\Psi^{(d)}_{i}$, for $i=1,\ldots,n$ into the tensors $\mathcal{T}^{d}$ (like in the previous paragraph), such that $\mathcal{T}^{d}_{i,:,\ldots,:}=\Psi^{(d)}_{i}$, for all $i=1,\ldots,n$. We now have for each degree a coupled partially symmetric CP decomposition as

[TABLE]

where the $\mathbf{c}_{i}$, for $i=1,\ldots,d$, are the $i$-th degree coefficients for each of the $r$ nonlinear functions $g_{k}$.

Remark that these coefficients were not required in the previous paragraphs when homogeneous polynomials were considered: in such cases the nonlinear functions $g_{k}$ are of the form $c_{k}t^{d}$, i.e., they differ only by a scaling factor, which can be assumed to be fully absorbed by $\mathbf{W}$. Also remark that there are redundancies in the representation (30): for example, an equivalent problem can be obtained if one rescales a coefficient vector $\mathbf{c}_{\delta}$ to a vector containing ones, in which case a rescaling has to take place on the remaining coefficients as well as on $\mathbf{W}$. Finally we want to mention that the framework of structured data fusion [43, 27] allows for computing tensor decompositions as in (30), where several tensors (and possibly matrices) are jointly decomposed while sharing factors, possibly while imposing structure on the factors.

Example 9**.**

Let us continue with Examples 1, 6 and 7. We have that $m=n=2$ and $r=3$, which does not guarantee a unique CP decomposition of $\mathcal{J}$ (under assumptions of genericity, see [24]). Indeed, if we compute a numerical CP decomposition of tensor $\mathcal{J}$, we find that, up to a relative norm-wise error $2.3546\times 10^{-16}$, $\mathcal{J}$ admits a CP decomposition with factors

[TABLE]

the columns of which are not scaled and permuted versions of the columns of $\mathbf{W}$, $\mathbf{V}$, $\mathbf{H}$.

It can be shown that the ‘structured CP’ approaches are able to correctly return the underlying factors $\mathbf{W}$, $\mathbf{V}$ and $\mathbf{H}$ (up to scaling and permutation invariances). For instance, the structured data fusion framework [43, 27] is able to compute the coupled and partially symmetric decomposition (30). This returns

[TABLE]

as well as computed values for the coefficient vectors of $g_{i}(x_{i})$, which are omitted here. It can be verified that $\tilde{\mathbf{W}}$ and $\tilde{\mathbf{V}}$ are scaled and permuted versions of $\mathbf{W}$ and $\mathbf{V}$.

Remark that if one uses $m=n=r=2$, both the structured and non-structured CP decomposition return the same decomposition (up to scaling and permutation of the columns of the factors). Indeed, in this case, uniqueness is guaranteed (generically), ensuring that the underlying factors are identifiable. This could be checked easily by generating a variation of the equations that we are decoupling where the third columns of $\mathbf{V}$ and $\mathbf{W}$ are removed, so that $g_{3}(x_{3})$ is not considered.

6.3 Linking $\mathcal{Q}$ and $\mathcal{T}^{1},\ldots,\mathcal{T}^{d}$ tensors

In this section, we show how $\mathcal{Q}$ and its CPD is connected with the tensors $\mathcal{T}^{s}$ and their joint decomposition (30). Let $(1,2)$-reshapings of the tensors $\mathcal{T}^{s}$ to be the third order tensors $\mathcal{T}^{s}_{(1,2)}\in\mathbb{R}^{n\times m\times m^{s-1}}$ defined as

[TABLE]

Then it is easy to see that the tensor $\mathcal{T}^{s}_{(1,2)}$ can be split into slices as shown in Fig. 5, where the $\Psi^{(s,j)}$ is the symmetric tensor corresponding to the $s$-th degree homogeneous part of the polynomial $f_{k}$.

By taking into account the definition (11) of the slices of the tensor $\mathcal{Q}$, we can easily see that $\mathcal{Q}$ can be constructed by stacking the tensors is equivalent to reshaping the tensors $\mathcal{T}^{s}_{(1,2)}$ along the third mode together, as shown in Fig. 6.

Remark that in Lemma 1 we see that the structure appearing in the CP decomposition of $\mathcal{Q}$ is closely connected to the simultaneous decomposition described in Section 6.2. Indeed, Lemma 1 can be alternatively deduced from (30), because the outer products of vectors become Kronecker products after reshaping.

7 Conclusions and perspectives

We have established a link between two tensorization approaches for decoupling multivariate polynomials [26, 24]: the tensor of Jacobian matrices [24] can be obtained by multiplying the coefficient-based tensor [26] by a Vandermonde-like matrix. As revealed by this connection, the two approaches have similar fundamental properties, such as equal tensor rank and uniqueness of the CP decomposition under conditions on the number and location of the sampling points.

The decoupling problem, however, is not equivalent to the CP decomposition of one of the tensors. This may lead to loss of uniqueness and identifiability of the CP decomposition, in the cases when the original decomposition is still unique. We have shown that by adding structure to the CP decomposition we can obtain equivalence between tensor decomposition and decoupling problems for polynomials. The structure can be imposed either as a joint decomposition of partially symmetric tensors, or can be imposed on rank-one factors. Numerical experiments confirm that using structured decompositions can restore uniqueness of the polynomial decoupling.

Although our results show that different tensor-based approaches are very closely related, let us make some remarks on applicability of the approaches and some future directions. For (differentiable) non-polynomial functions, the approach based on Jacobian matrices would be more appropriate, as it only uses evaluations of the derivatives of the functions. Coefficient-based approach seems more relevant in the case when the region of interest is unclear, or when some of the coefficients are missing or unreliable. In both cases, an interesting open question remains how to impose the structure directly on the rank-one components, without resorting to coupled tensor factorizations. Another important question is how to address the approximate decoupling problem, i.e., when we are dealing with noise (see [44] for results on the unstructured case).

8 Acknowledgements

This work was supported in part by the ERC AdG-2013-320594 grant “DECODA”, by the Flemish Government (Methusalem), by the ERC Advanced Grant SNLSID under contract 320378, and by the Fonds Wetenschappelijk Onderzoek – Vlaanderen under EOS Project no 30468160 and FWO projects G.0280.15N and G.0901.17N.