Beyond the Signs: Nonparametric Tensor Completion
via Sign Series

The signal plus noise model (3) is popular in tensor literature. Existing methods estimate the signal tensor based on low-rankness of $\Theta$ (Jain and Oh,, 2014; Montanari and Sun,, 2018). Common low-rank models include Canonical Polyadic (CP) tensors (Hitchcock,, 1927), Tucker tensors (De Lathauwer et al.,, 2000), and block tensors (Wang and Zeng,, 2019). While these methods have shown great success in signal recovery, tensors in applications often violate the low-rankness. Here we provide two examples to illustrate the limitation of classical models.

The first example reveals the sensitivity of tensor rank to order-preserving transformations. Let $\mathcal{Z}\in\mathbb{R}^{30\times 30\times 30}$ be an order-3 tensor with CP $\text{rank}(\mathcal{Z})=3$ (formal definition is deferred to end of this section). Suppose a monotonic transformation $f(z)=(1+\exp(-cz))^{-1}$ is applied to $\mathcal{Z}$ entrywise, and we let the signal $\Theta$ in model (1) be the tensor after transformation. Figure 1a plots the numerical rank (see Section 7.1) of $\Theta$ versus $c$. As we see, the rank increases rapidly with $c$, rending traditional low-rank tensor methods ineffective in the presence of mild order-preserving nonlinearities. In digital processing (Ghadermarzy et al.,, 2018) and genomics analysis (Hore et al.,, 2016), the tensor of interest often undergoes unknown transformation prior to measurements. The sensitivity to transformation makes the low-rank model less desirable in practice.

The second example demonstrates the inadequacy of classical low-rankness in representing special structures. Here we consider the signal tensor of the form $\Theta=\log(1+\mathcal{Z})$, where $\mathcal{Z}\in\mathbb{R}^{d\times d\times d}$ is an order-3 tensor with entries $\mathcal{Z}(i,j,k)={1\over d}\max(i,j,k)$ for $i,j,k\in\{1,\ldots,d\}$. The matrix analogy of $\Theta$ was studied by Chan and Airoldi, (2014) in graphon analysis. In this case neither $\Theta$ nor $\mathcal{Z}$ is low-rank; in fact, the rank is no smaller than the dimension $d$ as illustrated in Figure 1b. Again, classical low-rank models fail to address this type of tensor structure.

In the above and many other examples, the signal tensors $\Theta$ of interest have high rank. Classical low-rank models will miss these important structures. New methods that allow flexible tensor modeling have yet to be developed.

We develop a new model called sign representable tensors to address the aforementioned challenges. Figure 2 illustrates our main idea. Our approach is built on the sign series representation of the signal tensor, and we propose to estimate the sign tensors through a series of weighted classifications. In contrast to existing methods, our method is guaranteed to recover a wide range of low- and high-rank signals. We highlight two main contributions that set our work apart from earlier literature.

Statistically, the problem of high-rank tensor estimation is challenging. Existing estimation theory (Anandkumar et al.,, 2014; Montanari and Sun,, 2018; Cai et al.,, 2019) exclusively focuses on the regime of fixed $r$ growing $d$. However, such premise fails in high-rank tensors, where the rank may grow with, or even exceed, the dimension. A proper notion of nonparametric complexity is crucial. We show that, somewhat surprisingly, the sign tensor series not only preserves all information in the original signals, but also brings the benefits of flexibility and accuracy over classical low-rank models. The results fill the gap between parametric (low-rank) and nonparametric (high-rank) tensors, thereby greatly enriching the tensor model literature.

From computational perspective, optimizations regarding tensors are in general NP-hard. Fortunately, tensors sought in applications are specially-structured, for which a number of efficient algorithms are available (Ghadermarzy et al.,, 2018; Wang and Li,, 2020; Han et al.,, 2020). Our high-rank tensor estimate is provably reductable to a series of classifications, and its divide-and-conquer nature facilitates efficient computation. The ability to import and adapt existing tensor algorithms is one advantage of our method.

We also highlight the challenges associated with tensors compared to matrices. High-rank matrix estimation is recently studied under nonlinear models (Ganti et al.,, 2015) and subspace clustering (Ongie et al.,, 2017; Fan and Udell,, 2019). However, the problem for high-rank tensors is more challenging, because the tensor rank often exceeds the dimension when order $K\geq 3$ (Anandkumar et al.,, 2017). This is in sharp contrast to matrices. We show that, applying matrix methods to higher-order tensors results in suboptimal estimates. A full exploitation of the higher-order structure is needed; this is another challenge we address in this paper.

We use $\textup{sgn}(\cdot)\colon\mathbb{R}\to\{-1,1\}$ to denote the sign function, where $\textup{sgn}(y)=1$ if $y\geq 0$ and $-1$ otherwise. We allow univariate functions, such as $\textup{sgn}(\cdot)$ and general $f\colon\mathbb{R}\to\mathbb{R}$, to be applied to tensors in an element-wise manner. We denote $a_{n}\lesssim b_{n}$ if $\lim_{n\to\infty}a_{n}/b_{n}\leq c$ for some constant $c\geq 0$. We use the shorthand $[n]$ to denote the $n$-set $\{1,\ldots,n\}$ for $n\in\mathbb{N}_{+}$. Let $\Theta\in\mathbb{R}^{d_{1}\times\cdots\times d_{K}}$ denote an order-$K$ $(d_{1},\ldots,d_{K})$-dimensional tensor, and $\Theta(\omega)\in\mathbb{R}$ denote the tensor entry indexed by $\omega\in[d_{1}]\times\cdots\times[d_{K}]$. An event $E$ is said to occur “with very high probability” if $\mathbb{P}(E)$ tends to 1 faster than any polynomial of tensor dimension $d:=\min_{k}d_{k}\to\infty$. The CP decomposition (Hitchcock,, 1927) is defined by

where $\lambda_{1}\geq\cdots\geq\lambda_{r}>0$ are tensor singular values, $\bm{a}^{(k)}_{s}\in\mathbb{R}^{d_{k}}$ are norm-1 tensor singular vectors, and $\otimes$ denotes the outer product of vectors. The minimal $r\in\mathbb{N}_{+}$ for which (2) holds is called the tensor rank, denoted $\textup{rank}(\Theta)$.

Let $\Theta$ be the tensor of interest, and $\textup{sgn}(\Theta)$ the corresponding sign pattern. The sign patterns induce an equivalence relationship between tensors. Two tensors are called sign equivalent, denoted $\simeq$, if they have the same sign pattern.

The sign-rank is also called support rank (Cohn and Umans,, 2013), minimal rank (Alon et al.,, 2016), and nondeterministic rank (De Wolf,, 2003). Earlier work defines sign-rank for binary-valued tensors; we extend the notion to continuous-valued tensors. Note that the sign-rank concerns only the sign pattern but discards the magnitude information of $\Theta$. In particular, $\textup{srank}(\Theta)=\textup{srank}(\textup{sgn}\Theta)$.

Like most tensor problems (Hillar and Lim,, 2013), determining the sign-rank for a general tensor is NP hard (Alon et al.,, 2016). Fortunately, tensors arisen in applications often possess special structures that facilitate analysis. By definition, the sign-rank is upper bounded by the tensor rank. More generally, we have the following upper bounds.

Conversely, the sign-rank can be much smaller than the tensor rank, as we have shown in the examples of Section 1.

We provide several examples in Section 7.2, in which the tensor rank grows with dimension $d$ but the sign-rank remains a constant. The results highlight the advantages of using sign-rank in the high-dimensional tensor analysis. Propositions 1 and 2 together demonstrate the strict broadness of low sign-rank family over the usual low-rank family.

We now introduce a tensor family, which we coin as “sign representable tensors”, for the signal model in (3).

We show that the $r$-sign representable tensor family is a general model that incorporates most existing tensor models, including low-rank tensors, single index models, GLM models, and several hypergraphon models.

Accurate estimation of a sign representable tensor depends on the behavior of sign tensor series, $\textup{sgn}(\Theta-\pi)$. In this section, we show that sign tensors are completely characterized by weighted classification. The results bridge the algebraic and statistical properties of sign representable tensors.

For a given $\pi\in[-1,1]$, define a $\pi$-shifted data tensor $\bar{\mathcal{Y}}_{\Omega}$ with entries $\bar{\mathcal{Y}}(\omega)=(\mathcal{Y}(\omega)-\pi)$ for $\omega\in\Omega$. We propose a weighted classification objective function

where $\mathcal{Z}\in\mathbb{R}^{d_{1}\times\cdots\times d_{K}}$ is the decision variable to be optimized, $|\bar{\mathcal{Y}}(\omega)|$ is the entry-specific weight equal to the distance from the tensor entry to the target level $\pi$. The entry-specific weights incorporate the magnitude information into classification, where entries far away from the target level are penalized more heavily in the objective. In the special case of binary tensor $\mathcal{Y}\in\{-1,1\}^{d_{1}\times\cdots\times d_{K}}$ and target level $\pi=0$, the loss (4) reduces to usual classification loss.

Our proposed weighted classification function (4) is important for characterizing $\textup{sgn}(\Theta-\pi)$. Define the weighted classification risk

where the expectation is taken with respect to $\mathcal{Y}_{\Omega}$ under model (3) and the sampling distribution $\omega\sim\Pi$. Note that the form of $\textup{Risk}(\cdot)$ implicitly depends on $\pi$; we suppress $\pi$ when no confusion arises.

The results show that the sign tensor $\textup{sgn}(\Theta-\pi)$ optimizes the weighted classification risk. This fact suggests a practical procedure to estimate $\textup{sgn}(\Theta-\pi)$ via empirical risk optimization of $L(\mathcal{Z},\bar{\mathcal{Y}}_{\Omega})$. In order to establish the recovery guarantee, we shall address the uniqueness (up to sign equivalence) for the optimizer of $\textup{Risk}(\cdot)$. The local behavior of $\Theta$ around $\pi$ turns out to play a key role in the accuracy.

Some additional notation is needed. We use $\mathcal{N}=\{\pi\colon$ $\mathbb{P}_{\omega\sim\Pi}(\Theta(\omega)=\pi)\neq 0\}$ to denote the set of mass points of $\Theta$ under $\Pi$. Assume there exists a constant $C>0$, independent of tensor dimension, such that $|\mathcal{N}|\leq C$. Note that both $\Pi$ and $\Theta$ implicitly depend on the tensor dimension. Our assumptions are imposed to $\Pi=\Pi(d)$ and $\Theta=\Theta(d)$ in the high-dimensional regime uniformly where $d:=\min_{k}d_{k}\to\infty$.

The smoothness index $\alpha$ quantifies the intrinsic hardness of recovering $\textup{sgn}(\Theta-\pi)$ from $\textup{Risk}(\cdot)$. The value of $\alpha$ depends on both the sampling distribution $\omega\sim\Pi$ and the behavior of $\Theta(\omega)$. The recovery is easier at levels where points are less concentrated around $\pi$ with a large value of $\alpha>1$, or equivalently, when the cumulative distribution function (CDF) $G(\pi):=\mathbb{P}_{\omega\sim\Pi}[\Theta(\omega)\leq\pi]$ remains flat around $\pi$. A small value of $\alpha<1$ indicates the nonexistent (infinite) density at level $\pi$, or equivalently, when the $G(\pi)$ jumps at $\pi$. A typical case is $\alpha=1$ when the $G(\pi)$ has finite non-zero derivative in the vicinity of $\pi$. Table 1 illustrates the $G(\pi)$ for various models of $\Theta$ (see Section 5 Simulation for details).

We now reach the main theorem in this section. For two tensors $\Theta_{1},\Theta_{2}$, define the mean absolute error (MAE)

The result establishes the recovery stability of sign tensors $\textup{sgn}(\Theta-\pi)$ using optimization with population risk (5). The bound immediately shows the uniqueness of the optimizer for $\text{Risk}(\cdot)$ up to a zero-measure set under $\Pi$. We find that a higher value of $\alpha$ implies more stable recovery, as intuition would suggest. Similar results hold for optimization with sample risk (4) (see Section 4).

We conclude this section by applying Assumption 1 to the examples described in Section 3.1. For simplicity, suppose $\Pi$ is the uniform sampling for now. The tensor block model is $\infty$-globally smooth. This is because the set $\mathcal{N}$, which consists of distinct block means in $\Theta$, has finitely many elements. Furthermore, we have $\alpha=\infty$ for all $\pi\notin\mathcal{N}$, since the numerator in (7) is zero for all such $\pi$. The min/max hypergaphon model with a $r$-degree polynomial function is $1$-globally smooth because $\alpha=1$ for all $\pi$ in the function range except at most $(r-1)$ many stationary points.

Given a noisy incomplete tensor observation $\mathcal{Y}_{\Omega}$ from model (3), we cast the problem of estimating $\Theta$ into a series of weighted classifications. Specifically we propose the tensor estimate using the sign representation,

where $\hat{\mathcal{Z}}_{\pi}\in\mathbb{R}^{d_{1}\times\dots\times d_{K}}$ is the $\pi$-weighted classifier estimated at levels $\pi\in\mathcal{H}=\{-1,\ldots,-{1\over H},0,{1\over H},\ldots,1\}$,

Here $L(\cdot,\cdot)$ denotes the weighted classification objective defined in (4), where we have plugged $\bar{\mathcal{Y}}_{\Omega}=(\mathcal{Y}_{\Omega}-\pi)$ in the expression, and the rank constraint follows from Proposition 3. For the theory, we assume the true $r$ is known; in practice, $r$ could be chosen in a data adaptive fashion via cross-validation or elbow method (Hastie et al.,, 2009). Step (9) corresponds Figure 2c while (8) to Figure 2d .

The next theorem establishes the statistical convergence for the sign tensor estimate (9), which is an important ingredient for the final signal tensor estimate $\hat{\Theta}$ in (8).

Theorem 2 provides the error bound for the sign tensor estimation. Compared to the population results in Theorem 1, we here explicitly reveal the dependence of accuracy on the sample complexity and on the level $\pi$. The result demonstrates the polynomial decay of sign errors with $|\Omega|$. In particular, our sign estimate achieves consistent recovery using as few as $\tilde{O}(d_{\max}r)$ noisy entries.

Combining the sign representability of the signal tensor and the sign estimation accuracy, we obtain the main results on our nonparametric tensor estimation method.

Theorem 3 demonstrates the convergence rate of our tensor estimation. The bound (11) reveals three sources of errors: the estimation error for sign tensors, the bias from sign series representations, and the variance thereof. The resolution parameter $H$ controls the bias-variance tradeoff. We remark that the signal estimation error (12) is generally no better than the corresponding sign error (10). This is to be expected, since magnitude estimation is a harder problem than sign estimation.

In the special case of full observation with equal dimension $d_{k}=d,k\in[K]$, our signal estimate achieves convergence

Compared to earlier methods, our estimation accuracy applies to both low- and high-rank signal tensors. The rate depends on the sign complexity $\Theta\in\mathscr{P}_{\textup{sgn}}(r)$, and this $r$ is often much smaller than the usual tensor rank (see Section 3.1). Our result also reveals that the convergence becomes favorable as the order of data tensor increases.

We apply our method to the main examples in Section 3.1, and compare the results with existing literature. The numerical comparison is provided in Section 5.

The following sample complexity for nonparamtric tensor completion is a direct consequence of Theorem 3.

Our result improves earlier work (Yuan and Zhang,, 2016; Ghadermarzy et al.,, 2019; Lee and Wang,, 2020) by allowing both low- and high-rank signals. Interestingly, the sample requirements depend only on the sign complexity $r$ but not the nonparametric complexity $\alpha$. Note that $\tilde{\mathcal{O}}(d_{\max}r)$ roughly matches the degree of freedom of sign tensors, suggesting the optimality of our sample requirements.

This section addresses the practical implementation of our estimation (8) illustrated in Figure 2. Our sign representation of the signal estimate $\hat{\Theta}$ is an average of $2H+1$ sign tensors, which can be solved in a divide-and-conquer fashion. Briefly, we estimate the sign tensors $\mathcal{Z}_{\pi}$ (detailed in the next paragraph) for the series $\pi\in\mathcal{H}$ through parallel implementation, and then we aggregate the results to yield the output. The estimate enjoys low computational cost similar to a single sign tensor estimation.

For the sign tensor estimation (9), the problem reduces to binary tensor decomposition with a weighted classification loss. A number of algorithms have been developed for this problem (Ghadermarzy et al.,, 2018; Wang and Li,, 2020; Hong et al.,, 2020). We adopt similar ideas by tailoring the algorithms to our contexts. Following the common practice in classification, we replace the binary loss $\ell(z,y)=|\textup{sgn}z-\textup{sgn}y|$ with a surrogate loss $F(m)$ using a continuous function of margin $m:=z\textup{sgn}(y)$. Examples of large-margin loss are hinge loss $F(m)=(1-m)_{+}$, logistic loss $F(m)=\log(1+e^{-m})$, and nonconvex $\psi$-loss $F(m)=2\min(1,(1-m)_{+})$ with $m_{+}=\max(m,0)$. We implement the hinge loss and logistic loss in our algorithm, although our framework is applicable to general large-margin losses (Bartlett et al.,, 2006).

The rank constraints in the optimization (9) have been extensively studied in literature. Recent developments involve convex norm relaxation (Ghadermarzy et al.,, 2018) and nonconvex optimization (Wang and Li,, 2020; Han et al.,, 2020). Unlike matrices, computing the tensor convex norm is NP hard, so we choose (non-convex) alternating optimization due to its numerical efficiency. Briefly, we use the rank decomposition (2) of $\mathcal{Z}=\mathcal{Z}(\bm{A}_{1},\ldots,\bm{A}_{K})$ to optimize the unknown factor matrices $\bm{A}_{k}=[\bm{a}^{(k)}_{1},\ldots,\bm{a}^{(k)}_{r}]\in\mathbb{R}^{d_{k}\times r}$, where we choose to collect tensor singular values into $\bm{A}_{K}$. We numerically solve (8) by optimizing one factor $\bm{A}_{k}$ at a time while holding others fixed. Each suboptimization reduces to a convex optimization with a low-dimensional decision variable. Following common practice in tensor optimization (Anandkumar et al.,, 2014; Hong et al.,, 2020), we run the optimization from multiple initializations to locate a final estimate with the lowest objective value. The full procedure is described in Algorithm 1.

The brain dataset records the structural connectivity among 68 brain regions for 114 individuals along with their Intelligence Quotient (IQ) scores. We organize the connectivity data into an order-3 tensor, where entries encode the presence or absence of fiber connections between brain regions across individuals.

Figure 5 shows the MAE based on 5-fold cross-validations with $r=3,6,\ldots,15$ and $H=20$. We find that our method outperforms CPT in all combinations of ranks and missing rates. The achieved error reduction appears to be more profound as the missing rate increases. This trend highlights the applicability of our method in tensor completion tasks. In addition, our method exhibits a smaller standard error in cross-validation experiments as shown in Figure 5 and Table 2, demonstrating the stability over CPT. One possible reason is that that our estimate is guaranteed to be in $[0,1]$ (for binary tensor problem where $\mathcal{Y}\in\{0,1\}^{d_{1}\times\cdots d_{K}}$) whereas CPT estimation may fall outside the valid range $[0,1]$.

We next investigate the pattern in the estimated signal tensor. Figure 6a shows the identified top edges associated with IQ scores. Specifically, we first obtain a denoised tensor $\hat{\Theta}\in\mathbb{R}^{68\times 68\times 114}$ using our method with $r=10$ and $H=20$. Then, we perform a regression analysis of $\hat{\Theta}(i,j,\colon)\in\mathbb{R}^{144}$ against the normalized IQ score across the 144 individuals. The regression model is repeated for each edge $(i,j)\in[68]\times[68]$. We find that top edges represent the interhemispheric connections in the frontal lobes. The result is consistent with recent research on brain connectivity with intelligence (Li et al.,, 2009; Wang et al.,, 2017).

The NIPS dataset consists of word occurrence counts in papers published from 1987 to 2003. We focus on the top 100 authors, 200 most frequent words, and normalize each word count by log transformation with pseudo-count 1. The resulting dataset is an order-3 tensor with entry representing the log counts of words by authors across years.

Table 2 compares the prediction accuracy of different methods. We find that our method substantially outperforms the low-rank CP method for every configuration under consideration. Further increment of rank appears to have little effect on the performance. The comparison highlights the advantage of our method in achieving accuracy while maintaining low complexity. In addition, we also perform naive imputation where the missing values are predicted using the sample average. Both our method and CPT outperform the naive imputation, implying the necessity of incorporating tensor structure in the analysis.

We next examine the estimated signal tensor $\hat{\Theta}$ from our method. Figure 6b illustrates the results from NIPS data, where we plot the entries in $\hat{\Theta}$ corresponding to top authors and most-frequent words (after excluding generic words such as figure, results, etc). The identified pattern is consistent with the active topics in the NIPS publication. Among the top words are neural (marginal mean = 1.95), learning (1.48), and network (1.21), whereas top authors are T. Sejnowski (1.18), B. Scholkopf (1.17), M. Jordan (1.11), and G. Hinton (1.06). We also find strong heterogeneity among word occurrences across authors and years. For example, training and algorithm are popular words for B. Scholkopf and A. Smola in 1998-1999, whereas model occurs more often in M. Jordan and in 1996. The detected pattern and achieved accuracy demonstrate the applicability of our method.

In Section 1, we have provided a motivating example to show the sensitivity of tensor rank to monotonic transformations. Here, we describe the details of the example set-up.

The step 1 is to generate a rank-3 tensor $\mathcal{Z}$ based on the CP representation

where $\bm{a},\bm{b},\bm{c}\in\mathbb{R}^{30}$ are vectors consisting of $N(0,1)$ entries, and the shorthand $\bm{a}^{\otimes 3}=\bm{a}\otimes\bm{a}\otimes\bm{a}$ denotes the Kronecker power. We then apply $f(z)=(1+\exp(-cz))^{-1}$ to $\mathcal{Z}$ entrywise, and obtain a transformed tensor $\Theta=f(\mathcal{Z})$.

The step 2 is to determine the rank of $\Theta$. Unlike matrices, the exact rank determination for tensors is NP hard. Therefore, we choose to compute the numerical rank of $\Theta$ as an approximation. The numerical rank is determined as the minimal rank for which the relative approximation error is below $0.1$, i.e.,

We compute $\hat{r}(\Theta)$ by searching over $s\in\{1,\ldots,30^{2}\}$, where for each $s$, we (approximately) solve the least-square minimization using CP function in R package rTensor. We repeat steps 1-2 ten times, and plot the averaged numerical rank of $\Theta$ versus transformation level $c$ in Figure 1a.

In section 3.1, we have provided several tensor examples with high tensor rank but low sign-rank. This section provides more examples and their proofs. Unless otherwise specified, let $\Theta$ be an order-$K$ $(d,\ldots,d)$-dimensional tensor.

In fact, Example 6 is a special case of the following proposition.