Sparse and Low-rank Tensor Estimation via Cubic Sketchings

Introduction

The rapid advance in modern scientific technology gives rise to a wide range of high-dimensional tensor data . Accurate estimation and fast communication/processing of tensor-valued parameters are crucially important in practice. For example, a tensor-valued predictor which characterizes the association between brain diseases and scientific measurements becomes the point of interest . Another example is the tensor-valued image acquisition algorithm that can considerably reduce the number of required samples by exploiting the compressibility property of signals .

The following tensor estimation model is widely considered in recent literatures,

Here, $\mathscr{X}_{i}$ and $\epsilon_{i}$ are the measurement tensor and the noise, respectively. The goal is to estimate the unknown tensor $\mathscr{T}^{\ast}$ from measurements $\{y_{i},\mathscr{X}_{i}\}_{i=1}^{n}$. A number of specific settings with varying forms of $\mathscr{X}_{i}$ have been studied, e.g., tensor completion , tensor regression , multi-task learning , etc.

In this paper, we focus on the case that the measurement tensor can be written in a cubic sketching form. For example, $\mathscr{X}_{i}=\bm{x}_{i}\circ\bm{x}_{i}\circ\bm{x}_{i}$ or $\mathscr{X}_{i}=\bm{u}_{i}\circ\bm{v}_{i}\circ\bm{w}_{i}$, depending on whether $\mathscr{T}^{*}$ is symmetric or not. The cubic sketching form of $\mathscr{X}_{i}$ is motivated by a number of applications.

Interaction effect estimation: High-dimensional high-order interaction models have been considered under a variety of settings . By writing $\mathscr{X}_{i}=\bm{x}_{i}\circ\bm{x}_{i}\circ\bm{x}_{i}$, we find that the interaction model has an interesting tensor representation (see left panel of Figure 1) which allows us to estimate high-order interaction terms using tensor techniques. This is in contrast with the existing literature that mostly focused on pair-wise interactions due to the model complexity and computational difficulties. More detailed discussions will be provided in Section 5.

High-order imaging/video compression: High-order imaging/video compression is an important task in modern digital imaging with various applications (see right panel of Figure 1), such as hyper-spectral imaging analysis and facial imaging recognition . One could use Gaussian ensembles for compression such that each entry of $\mathscr{X}_{i}$ is i.i.d. randomly generated . In contrast, the non-symmetric cubic sketchings, i.e., $\mathscr{X}_{i}=\bm{u}_{i}\circ\bm{v}_{i}\circ\bm{w}_{i}$, reduce the memory storage from $O(np_{1}p_{2}p_{3})$ to $O(n(p_{1}+p_{2}+p_{3}))$ ($n$ is the sample size and $(p_{1},p_{2},p_{3})$ is the tensor dimension), but still preserve the optimal statistical rate. More detailed discussions will be provided in Section 6.

In practice, the total number of measurements $n$ is considerably smaller than the number of parameters in the unknown tensor $\mathscr{T}^{*}$, due to all kinds of restrictions such as time and storage. Fortunately, a variety of high-dimensional tensor data possess intrinsic structures, such as low-rankness and sparsity . This could highly reduce the effective dimension of the parameter and make the accurate estimation possible. Please refer to (3.2) and (6.2) for low-rankness and sparsity assumptions.

In this paper, we propose a computationally efficient non-convex optimization approach for sparse and low-rank tensor estimation via cubic-sketchings. Our procedure is two-stage:

obtain an initial estimate via the method of tensor moment (motivated by high-order Stein’s identity), and then apply sparse tensor decomposition to the initial estimate to output a warm start;

use a thresholded gradient descent to iteratively refine the warm start in each tensor mode until convergence.

Theoretically, we carefully characterize the optimization and statistical errors at each iteration step. The output estimate is shown to converge in a geometric rate to an estimation with minimax optimal rate in statistical error (in terms of tensor Frobenius norm). In particular, after a logarithmic number of iterations, whenever $n\gtrsim K^{2}(s\log(ep/s))^{\tfrac{3}{2}}$, the proposed estimator $\widehat{\mathscr{T}}$ achieves

with high probability, where $s$, $K$, $p$, and $\sigma^{2}$ are the sparsity, rank, dimension, and noise level, respectively. We further establish the matching minimax lower bound to show that (1.2) is indeed optimal over a large class of sparse low-rank tensors. Our optimality result can be further extended to the non-sparse case (such as tensor regression ) – to the best of our knowledge, this is the first statistical rate optimality result in both sparse and non-sparse low-rank tensor regressions.

The above theoretical analyses are non-trivial due to the non-convexity of the empirical risk function, and the need to develop some new high-order sub-Gaussian concentration inequalities. Specifically, the empirical risk function in consideration satisfies neither restricted strong convexity (RSC) condition nor sparse eigenvalue (SE) condition in general. Thus, many previous results, such as the one based on local optima analysis , are not directly applicable. Moreover, the structure of cubic-sketching tensor leads to high-order products of sub-Gaussian random variables. Thus, the matrix analysis based on Hoeffding-type or Bernstein-type concentration inequality will lead to sub-optimal statistical rate and sample complexity. This motivates us to develop new high-order concentration inequalities and sparse tensor-spectral-type bound, i.e., Lemmas 1 and 2 in Section 4.3. These new technical results are obtained based on the careful partial truncation of high-order products of sub-Gaussian random variables and the argument of bounded $\psi_{\alpha}$-norm , and may be of independent interest.

The literature on low-rank matrix estimation methods, e.g., the spectral method and nuclear norm minimization , is also related to this work. However, our cubic sketching model is by-no-means a simple extension from matrix estimation problems. In general, many related concepts or methods for matrix data, such as singular value decomposition, are problematic to apply in the tensor framework . It is also found that simple unfolding or matricizing of tensors may lead to suboptimal results due to the loss of structural information . Technically, the tensor nuclear norm is NP-hard to even approximate , and thus the method to handle tensor low-rankness is distinct from the matrix.

The rest of the paper is organized as follows. Section 2 provides preliminaries on notation and basic knowledge of tensor. A two-stage method for symmetric tensor estimation is proposed in Section 3, with the corresponding theoretical analysis given in Section 4. A concrete application to high-order interaction effect models is described in Section 5. The non-symmetric tensor estimation model is introduced and discussed in Section 6. Numerical analysis is provided in Section 7 to support the proposed procedure and theoretical results of this paper. Section 8 discusses extensions to higher-order tensors. The proofs of technical results are given in supplementary materials.

Preliminary

Throughout the paper, vector, matrix, and tensor are denoted by boldface lower-case letters (e.g., $\bm{x},\bm{y}$), boldface upper-case letters (e.g., $\bm{X},\bm{Y}$), and script letters (e.g., $\mathcal{X},\mathcal{Y}$), respectively. For any set $A$, let $|A|$ be the cardinality. The ${\rm diag}(\bm{x})$ is a diagonal matrix generated by $\bm{x}$. For two vectors $\bm{x}$ and $\bm{y}$, $\bm{x}\circ\bm{y}$ is the outer product. Define $\|\bm{x}\|_{q}:=(|x_{1}|^{q}+\cdots+|x_{p}|^{q})^{1/q}$. We also define the $l_{0}$ quasi-norm by $\|\bm{x}\|_{0}=\#\{j:x_{j}\neq 0\}$ and $l_{\infty}$ norm by $\max_{1\leq j\leq p}|x_{j}|$. Denote the set $\{1,2,\ldots,n\}$ by $[n]$. Let $\bm{e}_{j}$ be the canonical vectors, whose $j$-th entry equals to 1 and all other entries equal to zero. For any two sequences $\{a_{n}\}_{n=1}^{\infty},\{b_{n}\}_{n=1}^{\infty}$, we say $a_{n}=\mathcal{O}(b_{n})$ if there exists some positive constant $C_{0}$ and sufficiently large $n_{0}$ such that $|a_{n}|\leq C_{0}b_{n}$ for all $n\geq n_{0}$. We also write $a_{n}\asymp b_{n}$ if there exists $C,c>0$ such that $ca_{n}\leq b_{n}\leq Ca_{n}$ for all $n\geq 1$. Additionally, $C_{1},C_{2},\ldots,c_{1},c_{2},\ldots$ are generic constants, whose actual values may be different from line to line.

More generally, we may decompose a tensor as the sum of rank one tensors as follows,

where $\langle\mathcal{X},\mathcal{Y}\rangle=\sum_{i,j,k}\mathcal{X}_{ijk}\mathcal{Y}_{ijk}$. Clearly, $\|\mathcal{X}\|_{F}^{2}=\langle\mathcal{X},\mathcal{X}\rangle$. We also consider the following sparse tensor spectral norm,

By definition, $\|\mathcal{X}\|_{s}\leq\|\mathcal{X}\|_{op}$. Suppose $\mathcal{X}=\bm{x}_{1}\circ\bm{x}_{2}\circ\bm{x}_{3}$ and $\mathcal{Y}=\bm{y}_{1}\circ\bm{y}_{2}\circ\bm{y}_{3}$ are two rank-one tensors. Then it is easy to check that $\|\mathcal{X}\|_{F}=\|\bm{x}_{1}\|_{2}\|\bm{x}_{2}\|_{2}\|\bm{x}_{3}\|_{2}$ and $\langle\mathcal{X},\mathcal{Y}\rangle=(\bm{x}_{1}^{\top}\bm{y}_{1})(\bm{x}_{2}^{\top}\bm{y}_{2})(\bm{x}_{3}^{\top}\bm{y}_{3})$.

Symmetric Tensor Estimation via Cubic Sketchings

In this section, we focus on the estimation of sparse and low-rank symmetric tensors,

where $\bm{x}_{i}$ are random vectors with i.i.d. standard normal entries. As previously discussed, the tensor parameter $\mathscr{T}^{*}$ often satisfies certain low-dimensional structures in practice, among which the factor-wise sparsity and low-rankness commonly appear. We thus assume $\mathscr{T}^{*}$ is CP rank-$K$ for $K\ll p$ and the corresponding factors are sparse,

The CP low-rankness has been widely assumed in literature for its nice scalability and simple formulation . Different from the matrix factor analysis, we do not assume the tensor factors $\bm{\beta}_{k}^{*}$ here are orthogonal. On the other hand, since the low-rank tensor estimation is NP-hard in general , we will introduce an incoherence condition in the forthcoming Condition 3 to ensure that the correlation among different factors $\bm{\beta}_{k}^{*}$ is not too strong. Such a condition has been used in recent literature on tensor data analysis , compressed sensing , matrix decomposition , and dictionary learning .

Based on observations $\{y_{i},\mathscr{X}_{i}\}_{i=1}^{n}$, we propose to estimate $\mathscr{T}^{*}$ via minimizing the empirical squared loss since the close-form gradient provides computational convenience,

Clearly, (3.5) is a non-convex optimization problem. To solve it, we propose a two-stage method as described in the next two subsections.

Due to the non-convexity of (3.5), a straightforward implementation of many local search algorithms, such as gradient descent and alternating minimization, may easily get trapped into local optimums and result in sub-optimal statistical performance. Inspired by recent advances of spectral method (e.g., EM algorithm , phase retrieval , and tensor SVD ), we propose to evaluate an initial estimate $\{\eta_{k}^{(0)},\bm{\beta}_{k}^{(0)}\}$ via the method of moment and sparse tensor decomposition (a variant of high-order spectral method) in the following Steps 1 and 2, respectively. The pseudo-code is given in Algorithm 1.

Step 1: Unbiased Empirical Moment Estimator. Construct the empirical moment-based estimator ${\mathcal{T}}_{s}$,

Based on Lemma 4, ${\mathcal{T}}_{s}$ is an unbiased estimator of $\mathscr{T}^{\ast}$. The construction of (3.6) is motivated by the high-order Stein’s identity (; also see Theorem 7 for a complete statement). Intuitively speaking, based on the third-order score function of a Gaussian random vector $\bm{x}$: ${\mathcal{S}}_{3}(\bm{x})=\bm{x}\circ\bm{x}\circ\bm{x}-\sum_{j=1}^{p}(\bm{x}\circ\bm{e}_{j}\circ\bm{e}_{j}+\bm{e}_{j}\circ\bm{x}\circ\bm{e}_{j}+\bm{e}_{j}\circ\bm{e}_{j}\circ\bm{x})$, we can construct the unbiased estimator of $\mathscr{T}^{*}$ by properly choosing a continuously differentiable function in high-order Stein’s identity. See the proof of Lemma 4 for details.

Step 2: Sparse Tensor Decomposition. Based on the method of moment estimator obtained in Step 1, we further obtain good initialization for the factors $\{\eta_{k}^{(0)},\bm{\beta}_{k}^{(0)}\}$ via truncation and alternating rank-1 power iterations ,

Note that the tensor power iterations recover one rank-1 component per time. To identify all rank-1 components, we generate a large number of different initialization vectors, implement a clustering step, and choose the centroids as the estimates in the initialization stage. This scheme originally appears in tensor decomposition literature , although our problem setting and proof techniques are very different. This procedure is also very different from the matrix setting since the rank-1 component in singular value decomposition is mutually orthogonal, but we do not enforce the exact orthogonality here for $\mathscr{T}^{*}$.

where $\times_{2},\times_{3}$ are tensor multiplication operators defined in Section 2 and $T_{d}(\bm{x})$ is a truncation operator that sets all but the largest $d$ entries in absolute values to zero for any vector $\bm{x}$. It is noteworthy that the symmetry of ${\mathcal{T}}_{s}$ implies

This means the multiplications along different modes are the same. We run power iterations till its convergence, and denote $\bm{b}_{m}$ as the outcome. Finally, we apply $K$-means to partition $\{\bm{b}_{m}\}_{m=1}^{M}$ into $K$ clusters, let the centroids of the output clusters be $\{\bm{\beta}_{k}^{(0)}\}_{k=1}^{K}$, and calculate $\eta_{k}^{(0)}={\mathcal{T}}_{s}\times_{1}\bm{\beta}_{k}^{(0)}\times_{2}\bm{\beta}_{k}^{(0)}\times_{3}\bm{\beta}_{k}^{(0)}$ for $k\in[K]$.

2 Thresholded Gradient Descent

be the gradient function with respect to $\bm{B}$. Based on the detailed calculation in Lemma E.1, $\nabla_{\bm{B}}\mathcal{L}(\bm{B},\bm{\eta})$ can be written as

where $\{(\bm{B}^{\top}\bm{X})^{\top}\}^{3}$ and $\{(\bm{B}^{\top}\bm{X})^{\top}\}^{2}$ are entry-wise cubic and squared matrices of $(\bm{B}^{\top}\bm{X})^{\top}$. Define $\varphi_{h}(x)$ as the thresholding function with a level $h$ that satisfies the following minimal assumptions:

The initial estimates $\bm{\eta}^{(0)}$ and $\bm{B}^{(0)}$ will be updated by thresholded gradient descent in two steps summarized in Algorithm 2. It is noteworthy that only $\bm{B}$ is updated in Step 3, while $\bm{\eta}$ will be updated in Step 4 after finishing the update of $\bm{B}$.

Step 3: Updating $\bm{B}$ via Thresholded Gradient descent. We update $\bm{B}^{(t)}$ via thresholded gradient descent,

$\mu$ is the step size and $\phi=\sum_{i=1}^{n}y_{i}^{2}/n$ serves as an approximation for $(\sum_{k=1}^{K}\eta_{k}^{*})^{2}$ (see Lemma 15);

Step 4: Updating $\bm{\eta}$ via Normalization. We normalize each column of $\bm{B}^{(T)}$ and estimate the weight parameter as

The final estimator for $\mathscr{T}^{\ast}$ is

The evaluation of the gradient (3.7) requires $\mathcal{O}(npK^{2})$ operations at each iteration and can be computationally intense for large $n$ or $p$. To economize the computational cost, a stochastic version of thresholded gradient descent algorithm can be easily carried out by sampling a subset of summand functions (3.7) at each iteration. This will accelerate the procedure especially in the case of large-scale settings. See Section E.2 for details.

Theoretical Analysis

In this section, we establish the geometric convergence rate in optimization error and minimax optimal rate in statistical error of the proposed symmetric tensor estimator.

We first introduce the assumptions for theoretical analysis. Conditions 1-3 are on the true tensor parameter $\mathscr{T}^{*}$ and Conditions 4-5 are on the measurement scheme. Specifically, the first condition ensures the model identifiability for CP-decomposition.

The CP-decomposition in (3.2) is unique in the sense that if there exists another CP-decomposition $\mathscr{T}^{*}=\sum_{k=1}^{K^{\prime}}\eta_{k}^{*^{\prime}}\bm{\beta}_{k}^{*^{\prime}}\circ\bm{\beta}_{k}^{*^{\prime}}\circ\bm{\beta}_{k}^{*^{\prime}}$, it must have $K=K^{\prime}$ and be invariant up to a permutation of $\{1,\ldots,K\}$.

For technical purposes, we introduce the following conditions to regularize the CP-decomposition of $\mathscr{T}^{\ast}$. Similar assumptions were imposed in recent tensor literature, e.g., and Assumption 1.1 (A4) .

The CP-decomposition $\mathscr{T}^{*}=\sum_{k=1}^{K}\eta_{k}^{*}\bm{\beta}_{k}^{*}\circ\bm{\beta}_{k}^{*}\circ\bm{\beta}_{k}^{*}$ satisfies

for some absolute constants $C,C^{\prime}$, where $\eta_{\min}^{*}=\min_{k}\eta_{k}^{*}$ and $\eta_{\max}^{*}=\max_{k}\eta_{k}^{*}$. Recall that $s$ is the sparsity of $\bm{\beta}_{k}^{*}$.

In Condition 2, $R$ plays a similar role as a “condition number.” This assumption means that the tensor is “well-conditioned,” i.e., each rank-1 component is roughly of the same size.

As shown in the seminal work of , the estimation of low-rank tensors can be NP-hard in general. Hence, we impose the following incoherence condition.

The true tensor components are incoherent such that

where $R$ is the singular value ratio defined in (4.1) and $C^{{}^{\prime\prime}}$ is some small constant.

The preceding incoherence condition has been widely used in different scenarios in recent high-dimensional research, such as tensor decomposition , compressed sensing , matrix decomposition , and dictionary learning . It can be also viewed as a relaxation of orthogonality: if $\{\bm{\beta}_{1}^{*},\ldots,\bm{\beta}_{K}^{*}\}$ are mutually orthogonal, $\Gamma$ equals zero. We can show from both theory (Lemma 28 in the supplementary materials) and simulation (Section 7) that the low-rank tensor $\mathscr{T}^{*}$ induced by (3.2) satisfies the incoherence condition with high probability, if the component vectors $\bm{\beta}_{k}^{*}$ are randomly generated, say from Gaussian distribution.

We also introduce the following conditions on noise distribution.

The sample complexity condition is crucial for our algorithm especially in the initialization stage. Ignoring any polylog factors, Condition 5 is even weaker than the sparse matrix estimation case $(n\gtrsim s^{2})$ in .

2 Main Theoretical Results

Our main Theorem 1 shows that based on a proper initializer, the output of the proposed procedure can achieve optimal estimation error rate after a sufficient number of iterations. Here, we define the contraction parameter

and also denote $\mathcal{E}_{1}=4K\eta_{\max}^{*\tfrac{2}{3}}\varepsilon_{0}^{2}$ and $\mathcal{E}_{2}=C_{0}\eta_{\min}^{*-\tfrac{4}{3}}/16$ for some $C_{0}>0.$

Suppose Conditions 3-5 hold, $|\text{supp}(\bm{\beta}_{k}^{(0)})|\lesssim s$, and the initial estimator $\{\bm{\beta}_{k}^{(0)},\eta_{k}^{(0)}\}_{k=1}^{K}$ satisfy

with probability at least $1-\mathcal{O}(1/n)$. Assume the step size $\mu\leq\mu_{0}$, where $\mu_{0}$ is defined in (B.6). Then, the output of the thresholded gradient descent update in (3.9) satisfies:

For any $t=0,1,2,\ldots$, the factor-wise estimator satisfies

with probability at least $1-\mathcal{O}(tKs/n)$.

When the total number of iterations is no smaller than

there exists a constant $C_{1}$ (independent of $K,s,p,n,\sigma^{2}$) such that the final estimator $\widehat{\mathscr{T}}=\sum_{k=1}^{K}\eta_{k}^{(0)}\bm{\beta}_{k}^{(T^{*})}\circ\bm{\beta}_{k}^{(T^{*})}\circ\bm{\beta}_{k}^{(T^{*})}$ satisfies

with probability at least $1-\mathcal{O}(T^{*}Ks/n)$.

The error bound (4.3) can be decomposed into an optimization error $\mathcal{E}_{1}\kappa^{t}$ (which decays with a geometric rate as iterations) and a statistical error $\mathcal{E}_{2}\frac{\sigma^{2}s\log p}{n}$ (which does not decay as iterations). In the special case that $\sigma=0$, $\widehat{\mathscr{T}}$ exactly recover $\mathscr{T}^{\ast}$ with high probability.

The next theorem shows that Steps 1 and 2 of Algorithm 1 provides a good initializer required in Theorem 1.

Recall $\Gamma=\max_{1\leq k_{1}\neq k_{2}\leq K}|\langle\bm{\beta}_{k_{1}}^{*},\bm{\beta}_{k_{2}}^{*}\rangle|$. Suppose the number of initializations $L\geq K^{C_{3}\gamma^{-4}}$, where $\gamma$ is a constant defined in (B.3). Given that Conditions 1-4 hold, the initial estimator obtained from Steps 1-2 with a truncation level $s\leq d\leq Cs$ satisfies

Moreover, if the sample complexity condition 5 holds, then the above bound satisfies (4.2).

The upper bound of (4.6) consists of two terms that correspond to the approximation error of ${\mathcal{T}}_{s}$ to $\mathscr{T}^{\ast}$ and the incoherence among $\bm{\beta}_{k}^{\ast}$’s, respectively. Especially, the former converges to zero as $n$ grows while the latter does not.

The proof of Theorems 1 and 2 are postponed to Section B.1-B.2 in the supplementary materials. The combination of Theorems 1 and 2 immediately yields the following upper bound for the final estimator, which is one main result of this paper.

Suppose Conditions 1 – 5 hold, $s\leq d\leq Cs$. After $T^{*}$ iterations, there exists a constant $C_{1}$ not depending on $K,s,p,n,\sigma^{2}$, such that the proposed procedure yields

with probability at least $1-\mathcal{O}(T^{*}Ks/n)$, where $T^{*}$ is defined in (4.4).

The above upper bound turns out to match the minimax lower bound for a large class of sparse and low-rank tensors.

Consider the following class of sparse and low-rank tensors,

Suppose that $\{\mathscr{X}_{i}\}_{i=1}^{n}$ are i.i.d standard normal cubic sketchings with i.i.d. $N(0,\sigma^{2})$ noise in (3.1), $p\geq 20s$, and $s\geq 4$. We have the following lower bound result,

The proof of Theorem 4 is deferred to Section B.3 in the supplementary materials. Combining Theorems 3 and 4, we immediately obtain the following minimax-optimal rate for sparse and low-rank tensor estimation with cubic sketchings when $\log p\asymp\log(p/s)$:

The rate in (4.10) sheds light upon the effect of dimension $p$, noise level $\sigma^{2}$, sparsity $s$, sample size $n$ and rank $K$ to the estimation performance.

Recently, Li, Haupt, and Woodruff studied the optimal sketching for the low-rank tensor regression and gave an near-optimal sketching complexity with a sharp $(1+\varepsilon)$-worse-case error bound. Different from the framework of that focuses on a deterministic setting, we study a probabilistic model with random observation noises, propose a new algorithm, and studied the minimax optimal rate of estimation errors. In addition, considered different types of convex/non-convex algorithms for low-rank tensor regression with statistical assumptions. To our best knowledge, we are the first to achieve an optimal rate in estimation error based on polynomial-time algorithms for the tensor regression problem.

When the low-rank tensor $\mathscr{T}^{\ast}$ is not necessarily sparse, i.e.,

we can apply the proposed procedure with all the truncation/thresholding steps removed. If $n\geq\mathcal{O}(p^{3/2})$, we can use similar arguments of Theorems 1-3 to show that the estimator $\widehat{\mathscr{T}}^{\prime}$ satisfies

for any $\mathscr{T}^{*}\in\mathcal{F}_{p,K}$ with high probability. Furthermore, similar arguments of Theorem 4 imply that the rate in (4.11) is minimax optimal.

In addition, the analysis of gradient updates for the tensor case is significantly more complicated than the matrix case. First, it requires high-order concentration inequalities for the tensor case since the cubic-sketching tensor leads to high-order products of sub-Gaussian random variables (see Section 4.3 for details). The necessity of high-order expansions in the analysis of gradient updates for the tensor case also significantly increases the hardness of the problem. To ensure the geometric convergence, we need much more subtle analysis comparing to the ones in the matrix case .

3 Key Lemmas: High-order Concentration Inequalities

As mentioned earlier, one major challenge for theoretical analysis of cubic sketching is to handle heavy tails of high-order Gaussian moments. One can only handle up-to second moments of sub-Gaussian random variables by directly applying the Hoeffding’s or Bernstein’s concentration inequalities. Therefore, we need to develop the following high-order concentration inequalities as technical tools: Lemma 1 characterizes the tail bounds for the sum of sub-Gaussian products, and Lemma 2 provides the concentration inequalities for Gaussian cubic sketchings. The proofs of Lemmas 1 and 2 are given in Section A.2.

with probability at least $1-\delta$ for some constant $C$.

Note that in Lemma 1, each $\bm{X}_{i}$ does not necessarily have independent entries, even though $\{\bm{X}_{i}\}_{i=1}^{n}$ are independent matrices. Building on Lemma 1, Lemma 2 provides a generic spectral-type concentration inequality that can be used to quantify the approximation error of ${\mathcal{T}}_{s}$ introduced in Step 1 of the proposed procedure.

with probability at least $1-10/n^{3}-1/p$.

with probability at least $1-10/n^{3}-1/p$.

Here, $C$ is an absolute constant and $\|\cdot\|_{s}$ is the sparse tensor spectral norm defined in (2.3).

Application to High-order Interaction Effect Models

In this section, we study the high-order interaction effect model in the cubic sketching framework. Specifically, we consider the following three-way interaction model

Here $\bm{\xi}$, $\bm{\gamma}$, and $\bm{\eta}$ are coefficients for the main effect, pairwise interaction, and triple-wise interaction, respectively. More importantly, (5.1) can be reformulated into the following tensor form (also see the left panel of Figure 1)

We provide the following justification for assuming the tensorized coefficient $\mathcal{B}$ is low-rank and sparse. First, in modern applications, such as the biomedical research , the response is often driven by a small portion of coefficients and a small number of factors, leading to a highly entry-wise sparse and low-rank $\mathcal{B}$. Second, suggested that it is suitable to model entry-wise sparse and low-enough rank tensors as arising from sparse loadings. Therefore, we assume $\mathcal{B}$ is CP rank-$K$ with $s$-sparse factors:

where $K,s\ll p$. Then the number of parameters in (5.4), $K(p+1)$, is significantly smaller than $(p+1)^{3}$, the total number of parameters in the original three-way interaction effect model (5.1), which makes the consistent estimation of $\mathcal{B}$ possible in the high-dimensional case. In this case, (5.2) can be written as

By assuming $\bm{z}_{l}\overset{iid}{\sim}N_{p}(0,\bm{I}_{p})$, the high-order interaction effect model (5.2) reduces to the symmetric tensor estimation model (3.1), except one slight difference that the first coordinate of $\bm{x}_{l}$, i.e., the intercept, is always 1. To accommodate this difference, we only need to adjust the initial unbiased estimate in the above two-step procedure. Let

Then we construct the empirical moment-based initial tensor ${\mathcal{T}}_{s^{\prime}}$ as

For $i,j,k\neq 0$, ${\mathcal{T}}_{s^{\prime}[i,j,k]}={\mathcal{T}}_{s[i,j,k]}$, ${\mathcal{T}}_{s^{\prime}[i,j,0]}={\mathcal{T}}_{s[i,j,0]},{\mathcal{T}}_{s^{\prime}[0,j,k]}={\mathcal{T}}_{s[0,j,k]}$, and ${\mathcal{T}}_{s^{\prime}[i,0,k]}={\mathcal{T}}_{s[i,0,k]}$.

For $i\neq 0$, ${\mathcal{T}}_{s^{\prime}[0,0,i]}={\mathcal{T}}_{s^{\prime}[0,i,0]}={\mathcal{T}}_{s^{\prime}[i,0,0]}=\frac{1}{3}{\mathcal{T}}_{s[0,0,i]}-\frac{1}{6}(\sum_{k=1}^{p}{\mathcal{T}}_{s[k,k,i]}-(p+2)a_{i})$.

${\mathcal{T}}_{s^{\prime}}=\frac{1}{2p-2}(\sum_{k=1}^{p}{\mathcal{T}}_{s[0,k,k]}-(p+2){\mathcal{T}}_{s})$.

Lemma 5 shows that ${\mathcal{T}}_{s^{\prime}}$ is an unbiased estimator for $\mathcal{B}$.

The theoretical results in Section 4 imply the following upper and lower bounds for the three-way interaction effect estimation.

Suppose $\bm{z}_{1},\ldots,\bm{z}_{n}$ are i.i.d. standard Gaussian random vectors and $\mathcal{B}$ satisfies Conditions 1, 2 and 3. The output, denoted as $\widehat{\mathcal{B}}$, from the proposed Algorithms 1 and 2 based on ${\mathcal{T}}_{s^{\prime}}$ satisfies

with high probability. On the other hand, considering the following class of $\mathcal{B}$,

Non-symmetric Tensor Estimation Model

$\bm{B}_{1}=(\bm{\beta}_{11},\cdots,\bm{\beta}_{1K})$, $\bm{B}_{2}=(\bm{\beta}_{21},\cdots,\bm{\beta}_{2K})$, $\bm{B}_{3}=(\bm{\beta}_{31},\cdots,\bm{\beta}_{3K})$,

$\bm{U}=(\bm{u}_{1},\ldots,\bm{u}_{n}),\ \bm{V}=(\bm{v}_{1},\ldots,\bm{v}_{n}),\ \bm{W}=(\bm{w}_{1},\ldots,\bm{w}_{n})$, $\bm{\eta}=(\eta_{1},\ldots,\eta_{k})^{\top},\bm{y}=(y_{1},\ldots,y_{n})^{\top}.$

Then, the empirical risk function can be written compactly as

Since (6.3) is non-convex but fortunately tri-convex in terms of $\bm{B}_{1}$, $\bm{B}_{2}$, and $\bm{B}_{3}$, we develop a block-wise thresholded gradient descent algorithm as detailed below. The complete algorithm is deferred to Section 4 in the supplementary materials.

Construct the empirical moment-based estimator

to which sparse tensor decomposition is applied for initialization.

Lemma 17 shows that the gradient function for (6.3) with respect to $\bm{B}_{1}$ can be written as

where $\bm{D}=(\bm{B}_{1}^{\top}\bm{U})^{\top}*(\bm{B}_{2}^{\top}\bm{V})^{\top}*(\bm{B}_{3}^{\top}\bm{W})^{\top}\bm{\eta}-\bm{y}$ and $\bm{C}_{1}=(\bm{B}_{2}^{\top}\bm{V})^{\top}*(\bm{B}_{3}^{\top}\bm{W})^{\top}\odot\bm{\eta}^{\top}$. For $t=1,\ldots,T$, we fix $\bm{B}_{2}^{(t)},\bm{B}_{3}^{(t)}$ and update $\bm{B}_{1}^{(t+1)}$ via block-wise thresholded gradient descent,

where $\phi=\sum_{i=1}^{n}y_{i}^{2}/n$, $\mu$ is the step size, and $\bm{h}(\bm{B})=\sqrt{\frac{4\log np}{n^{2}}\{\bm{D}^{2}\}^{\top}\{\bm{C}^{2}\}}$. The updates of $\bm{B}_{2},\bm{B}_{3}$ are similar.

The theoretical analysis for the non-symmetric case is different from the symmetric one in two folds. First, the non-symmetric cubic sketching tensor is formed by three Gaussian vectors rather than one, which leads to many differences in the calculation of high-order moments. Second, the CP-decomposition of non-symmetric tensor $\mathscr{T}^{*}$ (6.2) forms a tri-convex optimization. At this point, the standard convex analysis for vanilla gradient descent could be applied given a proper initialization.

With the regularity conditions detailed in Section 4, we present the theoretical results for non-symmetric tensor estimation as follows.

Suppose Conditions 6 – 9 hold and $n\gtrsim(s\log(p_{0}/s))^{3/2}$, where $p_{0}=\max\{p_{1},p_{2},p_{3}\}$. For any $t=0,1,2,\ldots$, the output of Algorithm 4 satisfies

for some $0<\kappa<1$. When the total number of iterations is no smaller than $\log(\frac{n}{\sigma^{2}s\log p_{0}}\vee 1)/\log\kappa^{-1}$, the final estimator $\widehat{\mathscr{T}}$ satisfies

Consider the class of incoherent sparse and low-rank tensors $\mathcal{F}=\{\mathscr{T}:\mathscr{T}=\sum_{k=1}^{K}\bm{\beta}_{1k}\circ\bm{\beta}_{2k}\circ\bm{\beta}_{3k},\|\bm{\beta}_{i,k}\|_{0}\leq s\text{ for }i=1,2,3,k=1,\ldots,K\}$. If $\{\mathscr{X}_{i}\}_{i=1}^{n}$ are i.i.d standard normal cubic sketchings, $\epsilon\overset{iid}{\sim}N(0,\sigma^{2})$, $\min\{p_{1},p_{2},p_{3}\}\geq 20s$, and $s\geq 4$, we have

Theorems 5 and 6 imply that the proposed algorithm achieves a minimax-optimal rate of estimation error in the class of $\mathcal{F}$ as long as $\log(p_{0})\asymp\log(p_{0}/s)$.

Numerical Results

In this section, we investigate the effect of noise level, CP-rank, sample size, dimension, and sparsity on the estimation performance by simulation studies. We also investigate the numerical performance of the proposed algorithm when the incoherence assumption required in the theoretical analysis fails to hold.

In each setting, we generate $\mathscr{T}^{*}=\sum_{k=1}^{K}\bm{\beta}_{k}^{*}\circ\bm{\beta}_{k}^{*}\circ\bm{\beta}_{k}^{*}$, where $|\text{supp}(\bm{\beta}_{k}^{*})|=s$, the support of $\bm{\beta}_{k}^{*}$ is uniformly selected from $\{1,\ldots,p\}$, and the nonzero entries of $\bm{\beta}_{k}^{*}$ are drawn randomly from standard normal distribution. Then, we calculate $\eta_{k}^{*}\leftarrow\|\bm{\beta}_{k}^{*}\|_{2}^{3}$ and normalize $\bm{\beta}_{k}^{*}\leftarrow\bm{\beta}_{k}^{*}/\|\bm{\beta}_{k}^{*}\|_{2}$. The cubic sketchings $\{\mathscr{X}_{i}\}_{i=1}^{n}$ are generated as $\mathscr{X}_{i}=\bm{x}_{i}\circ\bm{x}_{i}\circ\bm{x}_{i}$ and $\bm{x}_{i}\overset{iid}{\sim}N(0,1)$. The noise satisfies $\{\epsilon_{i}\}_{i=1}^{n}\overset{iid}{\sim}N(0,\sigma^{2})$ or $\text{Laplace}(0,\sigma/\sqrt{2})$. Additionally, we adopt the following stopping rules in iterations: (1) the initialization iteration (Step 2 in Algorithm 1) is stopped if $\|\bm{b}_{m}^{(l+1)}-\bm{b}_{m}^{(l)}\|_{2}\leq 10^{-6}$; (2) the gradient update iteration (Step 3 in Algorithm 2) is stopped if $\|\bm{B}^{(T+1)}-\bm{B}^{(T)}\|_{F}\leq 10^{-6}$. The numerical results are based on 200 repetitions unless otherwise specified. The code was written in R and implemented on an Intel Xeon-E5 processor with 64 GB of RAM.

First, we consider the percentage of successful recovery in the noiseless case. Let $K=3$, $s/p=0.3$, $p=30$ or $50$, so that the total number of unknown parameters in $\mathscr{T}^{\ast}$ is $2.7\times 10^{4}$ or $1.25\times 10^{5}$. The sample size $n$ ranges from 500 to 6000. Each recovery is called “successful” if the relative error $\|\widehat{\mathscr{T}}-\mathscr{T}^{*}\|_{F}/\|\mathscr{T}^{*}\|_{F}<10^{-4}$. We report the average successful recovery rate in Figure 2.

We can see from Figure 2 that the empirical relation among successful recovery, dimension, and sample size is consistent with the theoretical results in Section 4.

We then move to the noisy case. Select $K=3$, $s/p=0.3$, $p\in\{30,50\}$, $\{\epsilon_{i}\}_{i=1}^{n}\overset{iid}{\sim}N(0,\sigma^{2})$. We consider two scenarios: (1) sample size $n$ = 6000, 8000, or 10000, $s/p=0.3$, the noise level $\sigma$ varies from 0 to 200; (2) noise level $\sigma=200$, sample size $n$ varies from 4000 to 10000, $p=30$, $s/p=0.1,0.3,0.5$. The estimation errors in terms of $\|\widehat{\mathscr{T}}-\mathscr{T}^{*}\|_{F}/\|\mathscr{T}^{*}\|_{F}$ in these two scenarios are plotted in Figures 4 and 4, respectively. These results show that the proposed procedure achieves a good performance – Algorithms 1 and 2 yield more accurate estimation with smaller variance $\sigma^{2}$ and/or large value of sample size $n$.

Next, we demonstrate that the low-rank tensor parameter $\mathscr{T}^{*}$ with randomly generated factors $\bm{\beta}_{k}^{*}$ satisfies the incoherence condition 3 with high probability. Set the CP-rank $K=3$ and the sparsity level $s/p=0.3$ with the dimension $p$ ranging from 10 to 2000. We compute the incoherence parameter $\Gamma$ defined in Condition 3. The left panel of Figure 5 shows that the incoherence parameter $\Gamma$ decays in a polynomial rate as $s$ grows, which matches the bound in Condition 3. Recall a theoretical justification on this point is also provided in Lemma 28.

We further examine the performance of the proposed algorithm when the incoherence condition required in the theoretical analysis fails to hold. Specifically, we set the CP-rank $K=3$, $p=30$, and the sparsity level $s/p=0.3$. We construct enormous copies of tensor parameter $\mathscr{T}^{*}_{j}$ with i.i.d. standard normal factor vectors $\bm{\beta}_{k}^{*}$. For each $\mathscr{T}^{*}_{j}$, we calculate the incoherence $\Gamma_{j}$ defined in Condition 3, then manually pick 40 $\mathscr{T}^{*}_{j^{\prime}}$ such that

In this way, we obtain a set of tensor parameters $\{\mathscr{T}^{*}_{j^{\prime}}\}$ with incoherence uniformly varying from 0 to 0.4. The right panel of Figure 5 plots the relative error for estimating $\mathscr{T}^{*}$ based on observations from cubic sketchings of $\mathscr{T}^{*}_{j^{\prime}}$ based on 1000 repetitions. We can see that the proposed algorithm achieves small relative errors even when the true factors are highly coherent.

Moreover, we consider a setting with Laplacian noise. Suppose $\{\epsilon_{i}\}_{i=1}^{n}\overset{iid}{\sim}Lap(\sigma)$ with density $f(x)=\frac{1}{\sigma}\exp(-2|x|/\sigma)$. With $n=3000$, $p=30$, and varying values of $\sigma$, the average estimation error and its comparison with Gaussian noise setting are provided in Figure 6. We note that the estimation errors under Laplace noise are slightly higher than those under Gaussian noise.

We also compare the estimation errors of initial and final estimators for different ranks and sample sizes. Set $K=3,p=30,s/p=0.3$ and consider the noiseless setting. It is clear from Figure 7 that the initialization error decays sufficiently, but does not converge to zero as sample size $n$ grows. This result matches our theoretical findings in Theorem 2: as discussed in Remark 5, the initial stage may yield an inconsistent estimator due to the incoherence among $\bm{\beta}_{k}$’s. We also evaluate and compare the estimation errors for both initial and final estimators. From the right panel of Figure 7, we can see that the final estimator is more stable and accurate compared to the initial one, which illustrates the merit of thresholded gradient descent step of the proposed procedure.

Finally, we compare the performance of the proposed method with the alternating least square (ALS)-based tensor regression method . We specifically consider two schemes for the initialization of ALS: (a) $\{\bm{\beta}_{k}^{(0)}\}$ are i.i.d. standard Gaussian (cold start), and (b) $\{\bm{\beta}_{k}^{(0)}\}$ are generated from the proposed Algorithm 1 (warm start). Setting $K=2$, $s/p=0.2$, $p=30$, $\{\epsilon_{i}\}_{i=1}^{n}\overset{iid}{\sim}N(0,200^{2})$, we apply both the proposed procedure and the ALS-based algorithm and record the average estimation errors with standard deviations for both initial and final estimators. From the result in Table 1, one can see the proposed algorithm significantly outperforms the ALS under both cold and warm start schemes. The main reason is pointed out in Remark 8: the cubic sketching setting possesses distinct aspects compared with the i.i.d. random Gaussian sketching setting, so that the method proposed by does not exactly fit here.

Discussions

This paper focuses on the third order tensor estimation via cubic sketchings. Moreover, all results can be extended to the higher-order case via high-order sketchings. To be specific, suppose

Then, one can similarly perform high-order sparse tensor decomposition and thresholded gradient descent to estimate $\mathscr{T}^{*}$. On the theoretical side, we can show if mild conditions hold and $n\geq C(\log n)^{m}(s\log p)^{m/2}$, the proposed procedure achieves

with high probability. The minimax optimality can be shown similarly.

References

Appendix A Proofs

We first introduce three lemmas to show that the empirical moment based tensors (3.6), (5.5), and (6.4) are all unbiased estimators for the target low-rank tensor in the corresponding scenarios. Detail proofs of three lemmas are postponed to Sections C.1.1, C.1.2 and C.1.3 in the supplementary materials.

For non-symmetric tensor estimation model (6.1) & (6.2), define the empirical moment-based tensor ${\mathcal{T}}$ by

Then ${\mathcal{T}}$ is an unbiased estimator for $\mathscr{T}^{*}$, i.e.,

The extension to the symmetric case is non-trivial due to the dependency among three identical sketching vectors. We borrow the idea of high-order Stein’s identity, which was originally proposed in . To fix the idea, we present only third order result for simplicity. The extension to higher-order is straightforward.

In general, the order-$m$ high-order score function is defined as

Interestingly, the high-order score function has a recursive differential representation

with ${\mathcal{S}}_{0}(\bm{x})=1$. This recursive form is helpful for constructing unbiased tensor estimator under symmetric cubic sketchings. Note that the first order score function ${\mathcal{S}}_{1}(\bm{x})=-\nabla\log p(\bm{x})$ is the same as score function in Lemma 26 (Stein’s lemma ). The proof of Theorem 7 relies on iteratively applying the recursion representation of score function (A.2) and the first-order Stein’s lemma (Lemma 26). We provide the detailed proof in Section B.4 for the sake of completeness.

In particular, if $\bm{x}$ follows a standard Gaussian vector, each order score function can be calculated based on (A.2) as follows,

Interestingly, if we let $G(\bm{x})=\sum_{k=1}^{K}\eta_{k}^{*}(\bm{x}^{\top}\bm{\beta}_{k}^{*})^{3}$, then

which is exactly $\mathscr{T}^{*}$. Connecting this fact with (A.1), we are able to construct the unbiased estimator in the following lemma through high-order Stein’s identity.

Consider the symmetric tensor estimation model (3.1) & (4.9). Define the empirical first-order moment $\bm{m}_{1}:=\frac{1}{n}\sum_{i=1}^{n}y_{i}\bm{x}_{i}$. If we further define an empirical third-order-moment-based tensor ${\mathcal{T}}_{s}$ by

Proof. Note that $y_{i}=G(\bm{x}_{i})+\epsilon_{i}$. Then we have

where ${\mathcal{S}}_{3}(\bm{x})$ is defined in (A.3). By using the conclusion in Theorem 7 and the fact (A.4), we obtain

since $\epsilon_{i}$ is independent of $\bm{x}_{i}$. This ends the proof. $\blacksquare$

Although the interaction effect model (5.1) is still based on symmetric sketchings, we need much more careful construction for the moment-based estimator, since the first coordinate of the sketching vector is always constant 1. We give such an estimator in the following lemma.

For interaction effect model (5.1), construct the empirical moment based tensor ${\mathcal{T}}_{s^{\prime}}$ as following

For $i,j,k\neq 0$, ${\mathcal{T}}_{s^{\prime}[i,j,k]}={\mathcal{T}}_{s[i,j,k]}$. And ${\mathcal{T}}_{s^{\prime}[i,j,0]}={\mathcal{T}}_{s[i,j,0]},{\mathcal{T}}_{s^{\prime}[0,j,k]}={\mathcal{T}}_{s[0,j,k]},{\mathcal{T}}_{s^{\prime}[i,0,k]}={\mathcal{T}}_{s[i,0,k]}$.

For $i\neq 0$, ${\mathcal{T}}_{s^{\prime}[0,0,i]}={\mathcal{T}}_{s^{\prime}[0,i,0]}={\mathcal{T}}_{s^{\prime}[i,0,0]}=\frac{1}{3}{\mathcal{T}}_{s[0,0,i]}-\frac{1}{6}(\sum_{k=1}^{p}{\mathcal{T}}_{s[k,k,i]}-(p+2)a_{i})$.

${\mathcal{T}}_{s^{\prime}}=\frac{1}{2p-2}(\sum_{k=1}^{p}{\mathcal{T}}_{s[0,k,k]}-(p+2){\mathcal{T}}_{s})$.

The ${\mathcal{T}}_{s^{\prime}}$ is an unbiased estimator for $\mathcal{B}$, i.e.,

A.2 Proofs of Lemmas 1 and 2: Concentration Inequalities

We aim to prove Lemmas 1 and 2 in this subsection. These two lemmas provide key concentration inequalities of the theoretical analysis for the main result. Before going into technical details, we introduce a quasi-norm called $\psi_{\alpha}$-norm.

The $\psi_{\alpha}$-norm of any random variable $X$ and $\alpha>0$ is defined as

Particularly, a random variable who has a bounded $\psi_{2}$-norm or bounded $\psi_{1}$-norm is called sub-Gaussian or sub-exponential random variable, respectively. Next lemma provides an upper bound for the $p$-th moment of sum of random variables with bounded $\psi_{\alpha}$-norm.

where $1/\alpha^{*}+1/\alpha=1$, $C_{1}(\alpha),C_{2}(\alpha)$ are some absolute constants only depending on $\alpha$.

If $0<\alpha<1$, (A.7) is a combination of Theorem 6.2 in and the fact that the $p$-th moment of a Weibull variable with parameter $\alpha$ is of order $p^{1/\alpha}$. If $\alpha\geq 1$, (A.7) follows from a combination of Corollaries 2.9 and 2.10 in . Continuing with standard symmetrization arguments, we reach the conclusion for general random variables. When $\alpha=1$ or 2, (A.7) coincides with standard moment bounds for a sum of sub-Gaussian and sub-exponential random variables in . The detailed proof of Lemma 6 is postponed to Section C.2.

When $0<\alpha<1$, by Chebyshev’s inequality, one can obtain the following exponential tail bound for the sum of random variables with bounded $\psi_{\alpha}$-norm. This lemma generalizes the Hoeffding-type concentration inequality for sub-Gaussian random variables (see, e.g. Proposition 5.10 in ), and Bernstein-type concentration inequality for sub-exponential random variables (see, e.g. Proposition 5.16 in ).

Proof. For any $t>0$, by Markov’s inequality,

where the last inequality is from Lemma 6. We set $t$ such that $\exp(-p)=C(\alpha)^{p}b^{p}(\sqrt{p}\|\bm{a}\|_{2}+p^{1/\alpha}\|\bm{a}\|_{\infty})^{p}/t^{p}$. Then for $p\geq 2$,

holds with probability at least $1-\exp(-p)$. Letting $\delta=\exp(-p)$, we have that for any $0<\delta<1/e^{2}$,

holds with probability at least $1-\delta$. This ends the proof. $\blacksquare$

The next lemma provides an upper bound for the product of random variables in $\psi_{\alpha}$-norm.

Suppose $X_{1},\ldots,X_{m}$ are $m$ random variables (not necessarily independent) with $\psi_{\alpha}$-norm bounded by $\|X_{j}\|_{\psi_{\alpha}}\leq K_{j}$. Then the $\psi_{\alpha/m}$-norm of $\prod_{j=1}^{m}X_{j}$ is bounded as

Proof. For any $\{x_{j}\}_{j=1}^{m}$ and $\alpha>0$, by using the inequality of arithmetic and geometric means we have

Since exponential function is a monotone increasing function, it shows that

From the definition of $\psi_{\alpha}$-norm, for $j=1,2,\ldots,m$, each individual $X_{j}$ has

Putting (A.8) and (A.9) together, we obtain

Therefore, we conclude that the $\psi_{\alpha/m}$-norm of $\prod_{j=1}^{m}X_{j}$ is bounded by $\prod_{j=1}^{m}K_{j}$. $\blacksquare$

Proof of Lemma 1. Note that for any $j=1,2,\ldots,m$, the $\psi_{2}$-norm of $\bm{X}_{j}^{\top}\bm{\beta}_{j}$ is bounded by $\|\bm{\beta}_{j}\|_{2}$ . According to Lemma 8, the $\psi_{2/m}$-norm of $\prod_{j=1}^{m}(\bm{X}_{j}^{\top}\bm{\beta}_{j})$ is bounded by $\prod_{j=1}^{m}\|\bm{\beta}_{j}\|_{2}$. Directly applying Lemma 7, we reach the conclusion. $\blacksquare$

Proof of Lemma 2. We first focus on the non-symmetric version and the proof follows three steps:

Truncate the first coordinate of $\bm{x}_{1i},\bm{x}_{2i},\bm{x}_{3i}$ by a carefully chosen truncation level;

Utilize the high-order concentration inequality in Lemma 20 at order three;

Show that the bias caused by truncation is negligible.

With slightly abuse of notations, we denote $a,x,y$ etc. as their first coordinate of $\bm{a},\bm{x},\bm{y}$ etc. Without loss of generality, we assume $p:=\max\{p_{1},p_{2},p_{3}\}$. By unitary invariance, we assume $\bm{\beta}_{1}=\bm{\beta}_{2}=\bm{\beta}_{3}=\bm{e}_{1}$, where $\bm{e}_{1}=(1,0,\ldots,0)^{\top}$. Then, it is equivalent to prove

Suppose $\bm{x}_{1}\sim\mathcal{N}(0,\bm{I}_{p_{1}}),\bm{x}_{2}\sim\mathcal{N}(0,\bm{I}_{p_{2}}),\bm{x}_{3}\sim\mathcal{N}(0,\bm{I}_{p_{3}})$ and $\{\bm{x}_{1i},\bm{x}_{2i},\bm{x}_{3i}\}_{i=1}^{n}$ are $n$ independent samples of $\{\bm{x}_{1},\bm{x}_{2},\bm{x}_{3}\}$. And define a bounded event $\mathcal{G}_{n}$ for the first coordinate and its corresponding population version,

we will prove that $M_{2}$ is negligible in terms of convergence rate of $M_{1}$.

Bounding $M_{1}$. For simplicity, we define $\bm{x}_{1}^{\prime}=\bm{x}_{1}|\mathcal{G},\ \bm{x}_{2}^{\prime}=\bm{x}_{2}|\mathcal{G},\ \bm{x}_{3}^{\prime}=\bm{x}_{3}|\mathcal{G}$, and $\{\bm{x}_{1i}^{\prime},\bm{x}_{2i}^{\prime},\bm{x}_{3i}^{\prime}\}_{i=1}^{n}$ are $n$ independent samples of $\{\bm{x}_{1}^{\prime},\bm{x}_{2}^{\prime},\bm{x}_{3}^{\prime}\}$. According to the law of total probability, we have

According to Lemma 22, the entry of $x_{1i}^{\prime}\bm{x}_{1i}^{\prime},x_{2i}^{\prime}\bm{x}_{2i}^{\prime},x_{3i}^{\prime}\bm{x}_{3i}^{\prime}$ are sub-Gaussian random variable with $\psi_{2}$-norm $M^{2}$. Applying Lemma 20, we obtain

where $\delta_{n,s}=((s\log(p/s))^{3}/n^{2})^{1/2}+(s\log(p/s)/n)^{1/2}$.

Putting the above bounds together, we obtain

By setting $M=2\sqrt{\log n/C_{2}}$, the bound of $M_{1}$ reduces to

Bounding $M_{2}$. From the definitions of $M_{2}$ and sparse spectral norm,

By the basic property of Gaussian random variable, we can show

where the last inequality holds for a large $M>0$. By the choice of $M=2\sqrt{\log n/C_{2}}$, we have $M_{2}\leq 208/C_{2}^{3/2}(\log n)^{\frac{3}{2}}/n^{2}$ for some constant $C_{2}$. When $n$ is large, this rate is negligible comparing with (A.10)

Bounding $M$: We put the upper bounds of $M_{1}$ and $M_{2}$ together. After some adjustments for absolute constant, it suffices to obtain

This concludes the proof of non-symmetric part. The proof of symmetric part remains similar and thus is omitted here. $\blacksquare$

Appendix B Additional Proofs for main results

Theorem 2 gives an approximation error upper bound for the sparse-tensor-decomposition-based initial estimator. In Step I of Section 3.1, the original problem can be reformatted to a version of tensor denoising:

The key difference between our model (B.1) and recent works is that $\mathcal{E}$ arises from empirical moment approximation, rather than the random observation noise considered in and . Next lemma gives an upper bound for the approximation error. The proof of Lemma 9 is deferred to Section C.3.

with probability at least $1-5/n$ for some uniform constant $C_{1}$.

Next we denote the following quantity for simplicity,

where $R$ is the singular value ratio, $K$ is the CP-rank, $s$ is the sparsity parameter, $\Gamma$ is the incoherence parameter and $C_{2}$ is uniform constant.

Next lemma provides theoretical guarantees for sparse tensor decomposition method.

Suppose that the symmetric tensor denoising model (B.1) satisfies Conditions 1, 2 and 3 (i.e., the identifiability, parameter space and incoherence). Assume the number of initializations $L\geq K^{C_{3}\gamma^{-4}}$ and the number of iterations $N\geq C_{4}\log\left(\gamma/\left(\frac{1}{\eta_{\min}^{*}}\|\mathcal{E}\|_{s+d}+\sqrt{K}\Gamma^{2}\right)\right)$ for constants $C_{3},C_{4}$, the truncation parameter $s\leq d\leq Cs$. Then the sparse-tensor-decomposition-based initialization satisfies

The proof of Lemma 10 essentially follows Theorem 3.9 in , we thus omit the detailed proof here. The upper bound in (B.4) contains two terms: $\frac{C_{4}}{\eta^{\ast}_{\min}}\|\mathcal{E}\|_{s+d}$ and $\sqrt{K}\Gamma^{2}$, which are due to the empirical moment approximation and the incoherence among different $\bm{\beta}_{k}$, respectively.

Although the sparse tensor decomposition is not optimal in statistical rate, it does offer a reasonable initial estimation provided enough samples. Equipped with (B.2) and Condition 2, the right side of (B.4) reduces to

with probability at least $1-5/n$. Denote $C_{0}=4\cdot 2160\cdot C_{1}C_{4}$. Using Conditions 3 and 5, we reach the conclusion that

with probability at least $1-5/n$. $\blacksquare$

B.2 Proof of Theorem 1: Gradient Update

Let $t\geq 0$ be an integer. Suppose Conditions 1-5 hold and $\{\bm{\beta}_{k}^{(t)},\eta_{k}\}$ satisfies the following upper bound

with probability at least $1-\mathcal{O}(K/n)$, where $\varepsilon_{0}=K^{-1}R^{-\tfrac{4}{3}}/2160$. As long as the step size $\mu$ satisfies

then $\{\bm{\beta}_{k}^{(t+1)}\}$ can be upper bounded as

with probability at least $1-\mathcal{O}(Ks/n)$.

In order to apply Lemma 11, we prove that the required condition (B.5) holds at every iteration step $t$ by induction. When $t=0$, by (4.2) and Condition 2,

holds with probability at least $1-\mathcal{O}(1/n)$. Since the initial estimator output by first stage is normalized, i.e., $\|\bm{\beta}_{k}^{(0)}\|_{2}=\|\bm{\beta}_{k}^{\ast}\|_{2}=1$, by triangle inequality we have

with probability at least $1-\mathcal{O}(1/n)$. Taking the summation over $k\in[K]$, we have

with probability at least $1-\mathcal{O}(K/n)$, which means (B.5) holds for $t=0$.

Suppose (B.5) holds at the iteration step $t-1$, which implies

for a sufficiently large $C_{0}$, we can obtain

By induction, (B.5) holds at each iteration step.

Now we are able to use Lemma 11 recursively to complete the proof. Repeatedly using Lemma 11, we have for $t=1,2,\ldots,$

with probability at least $1-\mathcal{O}(tKs/n)$. This concludes the first part of Theorem 1.

When the total number of iterations is no smaller than

the statistical error will dominate the whole error bound in the sense that

with probability at least $1-\mathcal{O}(T^{*}Ks/n)$.

The next lemma shows that the Frobenius norm distance between two tensors can be bounded by the distances between each factors in their CP decomposition. The proof of this lemma is provided in Section C.5.

Suppose $\mathscr{T}$ and $\mathscr{T}^{*}$ have CP-decomposition $\mathscr{T}=\sum_{k=1}^{K}\eta_{k}\bm{\beta}_{k}\circ\bm{\beta}_{k}\circ\bm{\beta}_{k}$ and $\mathscr{T}^{*}=\sum_{k=1}^{K}\eta_{k}^{*}\bm{\beta}_{k}^{*}\circ\bm{\beta}_{k}^{*}\circ\bm{\beta}_{k}^{*}$. If $|\eta_{k}-\eta_{k}^{*}|\leq c$, then

Denote $\widehat{\mathscr{T}}=\sum_{k=1}^{K}\eta_{k}\bm{\beta}_{k}^{(T^{*})}\circ\bm{\beta}_{k}^{(T^{*})}\circ\bm{\beta}_{k}^{(T^{*})}$. Combing (B.7) and Lemma 12, we have

with probability at least $1-\mathcal{O}(TKs/n)$. By setting $C_{1}=9C_{2}/4$, we complete the proof of Theorem 1. $\blacksquare$

B.3 Proofs of Theorems 4 and 6: Minimax Lower Bounds

Thus, if $\eta>0$, the moment generating function of $w^{(k,m_{1},m_{2})}-\frac{s}{2}$ satisfies

Here, (*) is due to $\eta>0$ and $\lfloor s/2\rfloor+1\geq s/2$. By setting $\eta=\log((p-s+1)/(8s))$, we have

Combining the two inequalities above, we have

Next we choose $M=\lfloor\exp(c_{0}/2\cdot sK\log(p/s))\rfloor$. Note that

which means there are positive probability that $\left\{\bm{\beta}^{(k,m)}\right\}_{\begin{subarray}{c}k=1,\ldots,K\\ m=1,\ldots,M\end{subarray}}$ satisfy

For the rest of the proof, we fix $\left\{\bm{\beta}^{(k,m)}\right\}_{\begin{subarray}{c}k=1,\ldots,K\\ m=1,\ldots,M\end{subarray}}$ to be the set of vectors satisfying (B.10).

where $\epsilon_{i}\overset{iid}{\sim}N(0,\sigma^{2}),\quad i=1,\ldots,n.$ Clearly, $\left(\bm{y}^{(m)},\bm{u},\bm{v},\bm{w}\right)$ follows a joint distribution, which may vary based on different values of $m$.

In this step, we analyze the Kullback-Leibler divergence between different distribution pairs:

Note that conditioning on fixed values of $\bm{u},\bm{v},\bm{w}$,

By the KL-divergence formula for Gaussian distribution,

Meanwhile, for any $1\leq m_{1}<m_{2}\leq M$,

By generalized Fano’s Lemma (see, e.g., ),

Finally we set $\lambda=\frac{c\sigma^{2}}{n}\log(p/s)$ for some small constant $c>0$, then

which has finished the proof of Theorem 6.

With $M=\exp(csK\log(p/s))$ and similar techniques as previous proof, one can show there exists positive possibility that

We then construct the following candidate symmetric tensors by blockwise design,

The rest of the proof essentially follows from the proof of Theorem 6. $\blacksquare$

B.4 Proof of Theorem 7: High-order Stein’s Lemma

The proof of this theorem follows from the one of Theorem 6 in . For the sake of completeness, we restate the detail here. Applying the recursion representation of score function (A.2), we have

Then, we apply the first-order Stein’s lemma (see Lemma 26) on function $G(\bm{x}){\mathcal{S}}_{2}(\bm{x})$ and obtain

Repeating the above argument two more times, we reach the conclusion. $\blacksquare$

Appendix C Proofs of Several Lemmas

In this subsection, we present the detail proofs of moment calculation, including non-symmetric case, symmetric case, and interaction model.

By the definition of $\{y_{i}\}$ in (6.1) & (6.2), we have

Each entry of $M$ can be calculated as follows

which implies $M=\bm{\beta}_{1}\circ\bm{\beta}_{2}\circ\bm{\beta}_{3}$. Combining with $n$ observations and $K$ components, we can obtain

C.1.2 Proof of Lemma 4

The last equation is due to $\|\bm{\beta}^{*}\|_{2}=1$.