High-Dimensional Tensor Regression With Oracle Properties
Authors: Wenbin Wang
Yu Shi
Ziping Zhao
Venue: ICML 2025 poster (2025)
Keywords: tensor regression
non-convex optimization
oracle properties
high-dimensional statistics
sparse
low rank
Primary area: general_machine_learning->evaluation

Abstract: The emergence of multi-dimensional data presents significant challenges for traditional regression models based on matrices or vectors, particularly in capturing multi-directional correlations. In response, tensor regression has been proposed as a powerful framework for modeling linear relationships among multi-dimensional variables. In this paper, we introduce a high-dimensional tensor-response tensor regression model under low-dimensional structural assumptions, such as sparsity and low-rankness. Assuming the underlying tensor lies within an unknown low-dimensional subspace, we consider a least squares estimation framework with non-convex penalties. Theoretically, we derive general risk bounds for the resulting estimators and demonstrate that they achieve the oracle statistical rates under mild technical conditions. To compute the proposed estimators efficiently, we introduce an accelerated proximal gradient algorithm demonstrating rapid convergence in practice. Extensive experiments on synthetic and real-world datasets validate the effectiveness of the proposed regression model and showcase the practical utility of the theoretical findings.

High-Dimensional Tensor Regression with Oracle Properties
Wenbin Wang'
Yu Shi
Ziping Zhao'
Abstract
in many domains (Han et al., 2022; Liu et al., 2022; Wang et al., 2024). For example, in image processing, tensor re- gression techniques have addressed critical tasks such as denoising (Zhang et al., 2021a), image inpainting (Bertalmio et al., 2001), and medical image analysis (Zhou et al., 2013; Li & Zhang, 2017). In recommendation systems, tensor regression leverages shared item information, substantially enhancing prediction accuracy and outperforming models that treat items independently (Zhang et al., 2021b). Addi- tionally, in spatio-temporal analysis, tensor regression has been developed to handle both forecasting and cokriging tasks (Bahadori et al., 2014; Yu & Liu, 2016; Yu et al., 2018; Su et al., 2020).
Tensor regression has emerged as powerful framework for modeling linear relationships among multi-dimensional variables by effectively capturing inherent cross-mode mode interactions within tensor-structured data. In this paper, we introduce high-dimensional tensor-response tensor regres- sion model under low-dimensional structural as- sumptions, such sparsity and low-rankness. Specifically, we assume that the underlying re- gression tensor lies within an unknown low- dimensional subspace and propose a general least squares estimation framework with non-conver penalties. Theoretically, we establish rigorous risk bounds for the resulting estimators, demonstrating that they attain the oracle statistical rates under mild technical conditions. To ensure computa- tional efficiency, we introduce a proximal gradient algorithm for solving the proposed non-convex op- timization problem. Extensive experiments con- datasets
Despite their wide applications. tensor regression models encounter significant challenges in estimation. This problem is more prominent in high-dimensional settings, where the number of parameters substantially exceeds the number of observations, which makes the model estimation ill-posed (Raskutti et al., 2019; Chen et al., 2019). To address such challenges, imposing structural assumptions that capture the underlying low-dimensional structures is critical, such as sparsity or low-rankness (Rabusseau & Kadri, 2016). Spar- sity refers to the phenomenon where most entries either exactly zero or near zero, which is commonly leveraged in fields such as recommendation systems (Lee, 2001) and compressed sensing (Chen et al., 2023). On the other hand, low-rankness indicates the rank is significantly smaller than its informative dimensions, which frequently applied in areas like image compression (Li & Li, 2010) and collab- orative filtering (Li et al., 2017). However, defining these low-dimensional structures for tensors is key challenge, as the concepts of sparsity and low-rankness have multiple nontrivial extensions in the tensor setting (Kolda & Bader, 2009). For instance, tensor sparsity can be defined either entry-wise or group-wise, such as at the fiber level (Raskutti et al., 2019) or the slice level (Zhang et al., 2019). Similarly, low-rankness can be imposed on different forms of tensors, such mode-wise (Chen et al., 2019) and slice-wise (Farias & Li, 2017; Luo & Zhang, 2024). Raskutti et al. (2019) in- vestigated all of the previously mentioned sparsity and low- rankness structures, establishing both general risk bounds and specific upper bounds in various scenarios. Alterna- tively, tensor decomposition techniques, such as Canonical Polyadic (CP) decomposition (Carroll & Chang, 1970) and
ducted on both synthetic and validate the effectiveness of the proposed regres- sion model and showease the practical utility of the theoretical findings.
1. Introduction
The tensor, a multi-diieensional array generalizing the ma- trix higher dimensions, has become useful tool in many data analysis areas (Kolda & Bader, 2009; Abraham et al., 2012; McConnell, 2014; McCullagh, 2018), including image analysis (Zhou et al., 2013; Li et al., 2018), biol- ogy (Hore et al., 2016), spectroscopy data (Amini et al., 2017), economics and finance (Li et al., 2015a; Wang et al., 2022), business (Hao et al., 2021: Bi et al., 2018), etc. Among tensor-based methods, the tensor regression is espe- cially useful for revealing linear relationships among high- dimensional variable sets and has been successfully applied
'School of Information Science and Technology, ShanghaiTech University, Shanghai 201210, China. Correspondence to: Ziping Zhao <zipingzhao shanghaitech.cdu..n.
Proceedings of the 4 International Conference on Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

High-Dimensionn Tensor Regression with Oracle Properties
Tucker decomposition (Tucker, 1966), offer another way to enforce structure. However, these methods often face chal- lenges due to their nonconvex nature of the optimization problems involved (Luo & Zhang, 2023). To further reduce the number of estimated parameters and improve model performance, several studies argue for decomposition estimators with regularizers (He et al., 2018; Ahmed et al., 2020: Xu, 2020).
matrices.
Tensor Unfolding A tensor A with order higher than 3 can be unfolded into an order-3 tensor along the (i, k)-th mode, denoted as A(j.)) E Rd xdlxxdxITessa d, This unfolding arranges the mode-(j, k) slices of A as the frontal slices of A(j,))- Specifically, its en- try satisfy [Ag.t)) 1.. = Ci where = + E L.= (1. - 1) II m/ dm- A tensor can also be unfolded into a matrix, which is also known as tensor matricization or flattening, The mode-(k) unfold- ing of A is denoted by A(a) E R d;, where each column corresponds to mode-(k) fiber of A. Specifi- cally, the entries of A(k) satisfy [A(K)] in.1 -fi..nde....d
In this paper, we investigate high-dimensional tensor-on- tensor regression model and propose to estimate the coeffi- cient with penalty regularizers. While convex methods, such as Lasso (Tibshirani, 1996) and nuclear norm minimization (Recht et al., 2010; Candes & Recht, 2012), are widely em- ployed due to their strong theoretical guarantees (Zhang et al., 2019; Raskutti et al., 2019), nonconvex approaches have recently garnered attention for their advantages in unbi- ased estimation and improved theoretical properties in high- dimensional settings. Although nonconvex regularizers are widely used, their benefits for high-dimensional tensor re- gression problems remain unclear. This paper closes this gap by proposing general and unifying estimation frame- work. In particular, we discuss class of decomposable non- convex penalty functions, including the smoothly clipped absolute deviation (SCAD) (Fan & Li, 2001) and the min- imax concave penalty (MCP) (Zhang, 2010). Leveraging these univariate penalty functions, we impose distinct low- dimensional structures spprrsity or low-rank constraints- on the regression coefficient tensor. To effectively solve the resulting optimization problems, we present a proximal gradient algorithm. Our analysis shows that these estimators enjoy oracle properties under mild assumptions which is faster than the estimators (Raskutti et al., 2019). Extensive numerical experiments on both synthetic data and real-world datasets validate the theoretical results and demonstrate the practical advantages and breadth of the proposed methods. Proofs are deferred to the Supplementary Materials.
where / = 1 + Σ syy (i. - (i.- I m. dm. Finally, tensor can also be reshaped into a vector, an operation commonly known as tensor vectorization. We denote the vectorization of A as vec (A), which is equivalent to vec- torizing its mode-1 unfolding: vec (A) = vec (A(1)-
We use calligraphic letters represent sets, such S. The support of set S is denoted by |8|. 1ls (-) denotes the projection onto the set S. For function J. f' denotes its derivative, V represents its gradient, and f denotes its Hessian. For functionals f (x) and g (x), we denote f (x) Z g (x) i 2 eg (x). I (x) S g (xx if f (x) S Cg and (x) = g (x) if og (x) S (a) S Cg (x) for some positive constants e and C.
3. Problem Formulation
In this section, we present a unified framework for high- dimensional tensor regression with nonconvex regularizers.
3.1. Tensor Regression
We consider the following generic tensor-on- on-ensso regres- sion model (Lock, 2018; Raskutti et al., 2019; Miao et al., 2022) with tensor coefficient AA' E IRd xdw:
2. Preliminary
Throughout the paper, we use boldface calligraphy letters for tensors, such as A. boldface uppercase letters for matri- ces, such as A, and standard lowercase letters for scalars, such as x. The order (or degree) of a tensor is defined as the number of modes it has; hence, matrices, vectors, and scalars are order-2, 1r and order-O tensors, respec- tively. For order-N tensor A E R4 x--xdy, where the mode-i dimension is d; i E [1....,N) position (i,,...,iN) iN) is denoted by dii...ix or [A] ii.....
y (X, E
where X E Rd xddxxxxdaar is the predictor variable with M S N. y xgdar+1--xud Rd 11 xxxxyyy is the noise. (5-) is the tensor contraction
the response variable, and &
product between two tensors, where its (iart1,... --th entry is defined as
N). the entry of A at
( A""")
di
das
I Iir...
0
Fibers and Slices By fixing all indices of a tensor A except for the k-th mode, we obtain mode-(A)) fibers, which are vectors. Mode-(j.A) k) slices are obtained by fixing all indices except for the j-th and k-th modes, resulting in
1M
Specifically, when M = N. the output is a scalar, in which case (1) becomes a collar-on-tensoon regression model (Zhou et al., 2013; Gui et al., 2016).
2

---

High- Highh-Dimensionl Tensor Regression with Oracle Properties
3.2. Proposed Problem
There exists a constant v > 0 such that the penally
function satisfies p', (t) 0 for all t 2 v;
Given collection of n samples { (Dvt), ()0))))00)
)}
(=1', which
9A (t) is symmetric, i.e., qx (-t) = 92 for all t;
is assumed to be generated from the observation model (I), our goal is to estimate the unknown coefficient tensor A' by solving the following regularized least squares estimation problem:
9 (t) is monotone and Lipschitz continuous, i.e. for
12 2 11, there exists nonnegative constant such that
(1(2))-9((1))) S O;
12-11
min AER oond xdx
1000) -
+ Rr (A)
9u (t) and q (t) pass through the origin, i.e. 9x(0) d (0) = 0;
2n 1=1
(2)
There exists positive constant X such that Ig (1) S A for all t.
where RA (.4) is a structural regularization term. For a
tensor A, ||||| F (A.A)² 2,
Such conditions are commonly employed in the analysis of nonconvex statistical estimation problems (Wang et al., 2014; Gui al., 2016; Fan al., 2018). The third condition introduces a curvature property that governs the degree of concavity of qx (-), and consequently, the level of noncon-
3.3. Nonconvex Regularization
In this paper, we consider class of regularizers Rx (A), which is defined based on nonconvex penalty function (-) with parameter À 2 0. Prototype examples of such regularizers include the SCAD (Fan & Li, 2001) and MCP (Zhang, 2010), and we introduce them in the following.
vexity of pA (-). For SCAD, these conditions are satisfied with V = ad and ( = a-I. while for MCP, we have v = aA and Ç 1.
The SCAD penalty function, proposed by (Fan & Li, 2001), is defined as
3.3.1. 1. SPARSITY REGULARIZATION
A/1,, 1
[o S X, for < |t| S ad;
A straightforward approach to induce sparsity within ten- sor A is to enforce entry-wise sparsity. This strategy draws inspiration from the well-known Lasso regression (Tibshi- rani, 1996) and has been extensively studied using convex regularization methods (Zhang et al., 2019; Raskutti et al., 2019). In contrast to prior works, we employ a nonconvex penalty. Specifically, the entry-wise sparsity regularizer is defined as
PA (t) =
2(a-1) +
(a+1)A²
for|t] > ad,
2
where > 2 is tuning parameter. The MCP penalty, introduced by (Zhang, 2010), is defined as
A|t| axa
12
PA (1)
2a
1>
2
di
dy Px (dt....
for some constant a Note that both the above non- convex functions Px (4) can be decomposed as P> (f) A |t| + qu (t), where |t] denotes the absolute value of t and
RA (A) Σ
In addition to promoting entry-wise sparsity, this regulariza- tion framework can be extended to incorporate other forms of sparsity in general N-mode tensor, such as fiber-wise sparsity and Slice-wise sparsity. Further formulations and discussions of these alternative tions are presented in Appendix A.
ạx (4) is concave function. Specifically, for SCAD, yx (1)
is given by
-(14 2)2 ax (t)
1A ad)
2(a-1) 1) 1)
regulariza-
+
2 2aa+1)x² -Md) 1 (1t) 2 ad),
3.3.2. -ANNRSSSSSS REGULARIZATION
where 1 (-) denotes the indicator function; for MCP, yx (1)
In addition to promoting sparsity, encouraging low-rank structure in tensors has demonstrated significant benefits in various applications (Nion & Sidiropoulos, 2010; Li & Li, 2010; Collins & Cohen, 2012; Semerci et al., 2014). There are multiple notions of rank for higher-order tensors (Kolda & Bader, 2009). In this section, we focus on mode-wise low-rankness, which involves penalizing the singular values of the mode- (k) unfoldings. However, the commonly used tensor nuclear norm penalty, which applies the & norm
is
t² 1 (1t) E aA) +
( aa)
(ai²- - \|tl| -Mtt
9x (t) = 2a
((t) 2 ad).
Furthermore, we impose the following regularity conditions on PACA (-) and ax (-), as detailed in Assumption 1.
Assumption 1. The functions px (t) and qx (t) satisfy the following conditions:
3

---

High- Hiigh-Dimensional Tensor Regression with Oracle Properties
to the singular values of the unfolded matrices, inevitably introduces non-negligible bias (Raskutti al., 2019). To alleviate this issue, we propose using nonconvex penalty applied to each singular value. Specifically, the mode-wise low-rankness regularizer is defined as
Assumptions 2 and 3, which characterize the curvature prop- erties of the empirical loss function L. analogous the classical RSC and RSM conditions commonly used in the literature on linear regression problems (Wang et al., 2014; Gui et al., 2016; Elenberg et al., 2018). If Assumptions 2 and 3 hold at the same time, it implies that L 2 μ. Lever- aging the methodology of (Candes & Tao, 2007), can be proven that the empirical loss function C satisfies both the RSC and RSM conditions with high probability.
min Sonddde.Tl, (Te d,3
RA (AA)
Σ
PA (o (A(k)),
where 0; (A())) denotes the i-th singular value of A(k): Additionally, this framework can be extended to accom- modate alternative forms of low-rankuess regularization. Further problem formulations and theoretical analyses of these extensions are provided in Appendix A.
Assumption 4. Assume that the concatenation of vectorized covariates from n samples, denoted as [vec(X ,..., ......ve(([2" voc(x("))T] follows multivariate Gaussian distribution with zero mean and covariance matrix Σ. We assume that there exist constants 1 such that the eigenvalues of Σ satisfy:
4. Main Theory
K S Tmin ( Omax
Smax
Σ) S K.
In this section, we present the theoretical results for the estimators from (2) under different scenarios and derive their corresponding estimation error bounds. We begin by presenting some preliminary assumptions.
Assumption 4 ensures that the covariance matrix is posi- tive definite and well-conditioned, which is crucial for avoid- ing degeneracies in the parameter space. Such conditions are commonly met in a range of statistical estimation prob- lems (Liu et al., 2014; Raskutti et al., 2019; Wei & Zhao, 2023). the covariates { (x())
4.1. Preliminaries
To facilitate the subsequent discussion, we define a local region as
are independent and identically distributed, the covariance matrix Σ block- diagonal, and Assumption 4 reduces to similar conditions the covariance matrix of each individual sample.
C E Rd xxxxxxxxx AD(
where D (-) denotes distance function. For example, in the discussion of entry-wise wise sparsity,
4.2. Statistical Error Analysis
D (B A')
B A lp:
In the following, we establish statistical error bounds for different regularizers Rx (-), providing insights into their performance under different conditions.
and in the discussion of mode-wise low- ranknesnssss we define
|1 (B) D D(B,A)) = || (B)
muc
4.2.1. SPARSITY REGULARIZATION
where FA* is a subspace associated with the unfolding of A' (to be defined explicitly later), F is its orthogonal complement, and || Ilaue denotes the nuclear norm, i.e.,
Before presenting a detailed analysis of the convergence rates associated with the entry-wise sparsity regularizer, we first introduce the notion of the oracle rate. The oracle rate refers to the statistical convergence rate achieved by the oracle estimator, which serves as an idealized benchmark under the assumption that the true parameter support is known a priori. This assumption allows the oracle estimator to attain the best possible theoretical performance.
||. | -.. (.())
Define the empirical loss function as L(A)
25 Σラ (0 (x 4> | We make the following assumptions on this loss function.
Assumption 2 (Restricted strong convexity (RSC)). For any A, 1 E C, there exists constant j satisfying p > such that
Assuming the true parameter support set for the entry-wise sparsity is S₁, the entry-wise sparse oracle estimator is de- fined
C (B) 2 C(A) + (VL(A), AA + 2
B
All
A 0
arg
min A.A.
L L(A).
Assumption 3 (Restricted smoothness (RSM)). For any
A.B there exists a constant L satisfying L such that
where S7 denotes the complement of the support set S ((1,..., af * 0). By the mean value theorem,
L
C(B) S C(A) + (VL A) B AIE.
it is easy to obtain that A satisfies 1|
A'Ilr S

---

High- Hiigh-Dimensional Tensor Regression with Oracle Properties
(AA)sss
S Now, we have the following
operator X (A): Rdixxxxy RuxdM-=Ixxxxd defined
result.
x(A) (x A)...., (x .Ai)
Theorem 5 (Entry-wise sparsity).
Suppose that Assump- do) and the
tions 1-4 hold. If >5 S. X x=
along with its adjoint operator
x² (E
true parameter tensor A' satisfies
R Ruxdarsa xaxxyy
Rd xuxxxxy
de- ① x00 where Ei=
which
fined as X" (E(1m))
min (ioody]e
(E S X) )) Aat 1Mr+1... Cu..... T TiM I
then the optimal solution A problem (2) satisfies
Then, we introduce the oracle statistical convergence rate of the mode-wise low-rank estimator, which is as- sumed to know the true singular value support S2 = { { oi [T> ((x)) * #0} in advance. Specifically, the mode-wise low-rank oracle estimator is defined as
A A 2
S V
Theorem 5 implies that the proposed estimator achieves the oracle rate under relatively mild assumptions. This perfor- mance is superior to that of the existing estimator, which uses an é norm penalty (Raskutti et al., 2019). In fact, condition (3) is referred as the minimum signal strength condition, and denotes the minimum signal strength. This condition is commonly employed in the analysis of noncon- vex penalized regression problems (Fan & Li, 2001; Zhang, 2010; Fan et al., 2018), and it is considered rather mild because in our analysis, we take v x À to be of the order log(ds drr))) which can be very small as the sample size increases.
A 0 arg
min L (AA..
By the mean value theorem, it is easy to obtain that A satisfies || 0 A' 2 [[VC( s.
v Sal Then,
have the following result.
Theorem 6 (Mode-wise low-rankness). Suppose that Assumptions /~4 hold. If 1 G,
X X
Visal
x (E(li))
where || llsp is the spectral (K) (k)
norm, and the true parameter tensor A satisfies
4.2.2. LOW-RANKNESS REGULARIZATION
2V S2
min (A())) 2 v
[[x* (3
Consider the matrix X E R of rank r. Its sin- gular value decomposition is given by X = VEV. where U mmmr contains the left singular vectors, V guxr contains the right singular vectors, and Σ Diag (a) (X),..., (X)) grxr is diagonal matrix of the singular values. We further define subspace F FX) and its orthogonal complement F= (X) as follows:!
nu
then the optimal solution A to problem (2) satisfies
T4
SS2
A A'
n
where Tk = 11: ([x (E(Im))
(8)) )1
F (X)==
W span (W) C U, span (WT)
C CV
F (X)
W span (W) 1 U, span (WT)& T 1
This result suggests that, with an appropriately chosen regularization parameter À and provided that the smallest nonzero singular value is sufficiently large, the estimator will achieve the same rate as the oracle estimator.
where span(W) (W) denotes the subspace spanned by the columns of W. The projection operators onto the sub- space F and its orthogonal complement Frd are defined as follows:
5. Optimization Algorithm
T (X) U XVVT,
In this section, we present an accelerated proximal gradient algorithm, adapted from (Yao al., 2017), for solving the estimation problem (2). The algorithm iteratively combines gradient descent step on the loss function L(A) with proximal step on the regularization term R (A). To accelerate convergence, we incorporate an extrapolation step based on Nesterov's s accelerated gradient method (Nesterov, 2013).
1 (X) (I--U'')x((VVV))
where I denotes the identity matrix with contextually ap- propriate dimensions. Additionally, we introduce the linear
For brevity, we adopt the shorthand notations F and F when the dependence on X clear from the context.
5

---

High- Dimennionn Tensor Regression with Oracle Properties
Algorithm 1 Accelerated Proximal Gradient Algorithm
we adopt q balance the convergence speed and bility.
Require: E (0, 2), 0 E (0, L), X;
Au = A = 0:
Optimality Conditions Recall that the non-convex penalty px(t) can be decomposed as the sum of convex term A/t] and a concave component yx (t). Hence, the prob- lem (2) can be reformulated as
2 1 1;
repeat
}/ A (42 (A: A:-1);
A: = max, max;=max(1, g-.).... C (Ai)) + Rr (Ai);
L (V:) +Ra (y:)= S 4 then
min L(A) + (x + A|IAII
V.- =V,;
else
where (x(A) is the concave component of the Rx (A), and ||||| is a generic convex norm. For instance, in the entry-wise sparsity regularizer case,
V1 As:
10
end
V. Ar + 142 (A, - Ai-1);
d2
dy
12
3,- = V.- - VVVV...
2x(A)
aa (ail...;h))
13
A P (21)
EN =1
14
until WA (At) S
and ||||| is the & norm. In the mode-wise low-rankness regularizer case,
Output: Ar+1
min(Ie.I1: jonk d,}
Accelerated Proximal Gradient Algorithm Define the proximal operator associated with R-)-) as
2x
Σ
(A(o)
1 | 2
xrox (V) = arg min
VI RA
and ||||| is the nuclear norm of mode-(k) unfolded matrix.
After this decomposition, the optimization task simplifies considerably if the term (AA is omitted, reducing the classical linear programming problem. Motivated by this observation, the objective function can be reformulated as
Let { be the sequence of iterates. At the iteration I, the algorithm updates x through the following steps:
y: = Xx, Os (X, x))1), -
X++1 = x (())))))))))))
min ((A) |||||
where 0, 1-1 and 17 is the stepsize (Beck& & Teboulle, r+2 2009; Nesterov, 2013). When 0, 0. the algorithm reduces the standard proximal gradient algorithm. However, in the non-convex tensor setting, the extrapolation step V, can pre- vent guaranteed sufficient decrease of the objective func- tion in each iteration. To address this issue, Yao et al. (2017) proposed to evaluate the objective value before the proximal step, instead of checking after it as was previously done, demonstrating that this modification guarantees a sufficient decrease under certain mild conditions. To effectively imple- ment the algorithm,
1-1
where Ž(AA = Qx can be served as surrogate function and ||||| as a new convex penalty.
Since A is the exact global solution to the optimization prob- lem (2). By the Karush- Karush-Kuhn-- satisfies the following first-order optimal condition:
(KKT) conditions, A
Vī(Ā) + A = -0,
(4)
where タ =B||A]] is subgradient of || Equivalently, for all A.
need choose an appropriate step size 17. The parameter 7 is related to the Lipschitz constant L of (-) (i.e., IIV C (x) - V (>))lr < L|| || - - pp).) Choosing < 1/L ensures that the gradient step does not overshoot, which is critical for the convergence of the algorithm. The complete procedure summarized in Algo- rithm In line 5 of Algorithm q is an adjustable param- eter that governs nonmonotonic behavior of the objective sequences. Setting q 0 represents the simplest approach enforcing strictly monotone objective sequence, while larger ! permits the objective function C (Q))) + RA (2r) to occasionally increase. This flexibility is inspired by the Barzilai- Borreii--Boeennne scheme for unconstrained smooth mini- mization (Grippo& Sciandrone, 2002). In our experiments,
A-A,Vū(A) + >²) < 5.
(5)
Based on the optimality condition in (5), we measure the suboptimality of A using
wx (A) min GE
VŌ(A) + +
where ||. denote the dual norm of ||-I|- We say A is e-optimal solution to (2) if (x(A) < é with small positive number. Intuitively, when A is the exact global optimum, (AA S 0 by the KKT condition (4); otherwise, if A is optimal then (x(A) remains small but slightly positive.
6

---

High- Bich--Dimessonal Tensor Regression with Oracle Properties
0.8
h-1000
0.7
n=2000
0.9
n=2000
077
n-5000
n-3000
0.8
0.5
0.4
0.6
0.4
0.2
0.5
0.3
0.1
0.2
1024
1538
3000
[5]
(b)
(c)
Figure 1: Entry-wise sparsity regularizer with the error bars of MSFE + standard deviation.
no: 1000
n=:
n-2000
n-2000
0.
d-12
P-0000
1.5
n-3000
0.8
d-16
d-20
0.8
MSPE
0.8
0.4
0.5
0.
0
0
0.2
3000
[5
(b)
(c)
Figure 2: Mode-wise lowrankness regularizer with the error bars of MSFE = standard deviation.
6. Numerical Experiments
aiming to minimize the estimation error. All the reported results are averaged on 100 Monte Carlo realizations to ensure statistical robustness.
In this section, evaluate the performance of the proposed tensor regression model with various regularization schemes, as well as the optimization algorithm. In all experiments, the SCAD penalty is employed as the nonconvex regular- izer. The estimation performance is measured by the Mean Squared Frobenius norm Error (MSFE) and the Root Mean Square Error (RMSE). Specifically, the MSFE is defined as
6.1. Synthetic Data
We first evaluate our proposed estimator through compre- hensive synthetic experiments. We generate a set of indepen- dent covariate tensors (x(() Ji=1+ each entry independently drawn from standard Gaussian distribution. The response variables are obtained according model (1), with additive independent Gaussian noise having zero mean and variance parameterized by n For all synthetic experiments, we em- ploy 3rd-order tensors A E ddxdxd For Figures 1 and 2, the noise parameter is set to 0.1.
MSFE
A"
A
1==
where A is the estimated results and A is the true value, and the RMSE is defined as
Figure shows the estimation performance when employing the entry-wise sparsity regularizer, varying the dimension d, the number of nonzero entries |S,] and the sample size respectively. In Figure la and 1b. three lines correspond to sample size {1000, 2000, 3000). The proportion of non-zero entries s* /d is set to 0.5 for Figures la and Ic. Figure 1a demonstrates that a large sample size n consis- tently lower estimation error. In contrast, Figure 1b shows that higher sparsity levels raise estimation error.
y0 (x l.
RMSE
The tuning parameter À and the hyperparameter within the SCAD penalty are selected via ten-fold cross-validation,
II the whole section, use A to denote generic tensor estimator, which be the estimation results obtained various algorithms estimators.

---

HlihhDimensiinn Dimensional Tensor Regression with Oracle Properties
(a) nonconvex method 3D results
(b) convex method 3D results
(c) contour the tensor slice
Noroged
14
(d) nonconvex method 2D results
(c) convex method 2D results
(f) averaged thresholded value counts
Figure 3: Visualization of the estimated tensors with nonconvex and convex methods.
Table Comparisons between proposed nonconvex regularizers and convex regularizers.
Siructures
Methods
Synthetic Data
Real world Data
MSFE
RMSE
MSFE
MPRE
Nonconvex
0.4042 + 0.0201
+ 0.0021
134.5864 + 11.2950
7.6072 1
entry-wise sparsity
16 16 x
2048
0.1
Convex
0.6936 + 0.0297
0.1004 0.0023
144.7100 14.9947
7.7496 0.0457
Nonconvex
0.5482 1t 0.0395
0.1002 1: 0.0012
35.5536 1: 1.4889
1,0330 1 0.0022
made wise low raannees
16 x
Convex
1.7411 +: 0.0953
0.1096 4 0.0020
41,2719 4: 3.5079
1.1027 4: 0.0024
Figure illustrates the performance of the mode-wise low- rankness penalty. In Figure 2b, the w-axis [Sa] represent the rank of the mode-(k) unfolded matrix. Figures 2a and 2c set the rank of at 5, while Figure 2b and 2c fix the dimension of each mode 16. Each figure displays three distinet lines that correspond to the estimation errors for sample sizes of (1000, 2000, 3000).
counts of values above or below the threshold.
Table I compares the performance of our proposed non- convex regularizers against traditional convex regularizers (Raskutti et al., 2019). For sparsity, we set the v) 0.1, and for low-rankness,
1. We configure the tensor di- mension such A E ddxdxx or A E ddddddd with d 16. = 20. In all settings, our proposed regularizers achieve lower MSFE and RMSE than convex methods, aligning with our theoretical analysis.
In Figure 3, we plot each point in the matrix using its indices as coordinates, with the corresponding value on the z-axis. threshold color points to blue their absolute error is below this threshold, and red indicates otherwise. Figure 3a and 3b compare nonconvex and convex methods the tensor slice based on the mode-wise low-rank structure, where the size is 10 x 10, the sample number n = 1000 and the noise parameter 77 = 0.1. Figure 3c shows the contour of the original tensor slice, with the estimation threshold set to 10 Figures 3d and 3e visualize the thresholded results of the estimation in 2D, and Figure 3f illustrates average
6.2. Real-world Datasets
We validate our method on ImageNet 2012 dataset (Rus- sakovsky et al., 2015) for image denoising using a 3rd- order tensor AA E [66466433 with n = 4000 samples. In addition to MSFE, we also report Mean Prediction Rela-
- ->"" ll
tive Error (MPRE), defined as E
where
||
8

---

High--Dimensionl Tensor Regression with Oracle Properties
(a) original
(b) noisy
(c) convex
(d) nonconvex (prop.)
Figure 4: The original image, noisy image, and denoised images using convex and nonconvex methods.
V" (A", and 5 is the i-th prediction. Figure 4 shows the estimated image, and Table presents the com- parative results.
care, environmental monitoring. and multimedia), we do identify any specific foreseeable consequences that require explicit attention this time.
7. Conclusions and Future Work
References
Abraham, R., Marsden, J. E.,, and Ratiu, T. Manifolds, lensor analysis, and applications, volume 75. Springer Science & Business Media, 2012.
In this paper, we propose a comprehensive framework for tensor regression estimation using nonconvex regularizers. Our findings demonstrate that estimators employing noncon- vex regularizers exhibit faster convergence rates compared to those with convex regularizers. Furthermore, we show that under several mild conditions, our proposed estima- tor possesses the oracle property. Extensive experimental results validate our theoretical claims, showcasing a close alignment between the theoretical predictions and the ob- served numerical performance of our estimators. Currently, we are limited to applying regularization regularizers to tensor regression models. It would be desirable to derive some theoretical guarantees for alternative methods that capture structure in the tensor regression models, such as tensor decomposition; this is the aim our future work. To conclude, our work effectively bridges the gap between practical applications and theoretical analysis of tensor-on- tensor regression with nonconvex regularizers. To the best of our knowledge, this is the first work to obtain the ora- cle statistical rate of convergence for the tensor regression problem.
Ahmed, T., Raja, H., and Bajwa, W. U. Tensor regression using low-rank and sparse tucker decompositions. SIAM Journal on Mathematics of Data Science, 2(4):944-96, 2020.
Amini, A. A., Levina, E., and Shedden, K. A. Struc- tured regression models for high-dimensional spatial spec- troscopy data. Electronic Journal of Statistics, 11:4151- 4178, 2017.
Bahadori, M. T., Yu, Q. R., and Liu, Y. Fast multivariate spatio-temporal analysis via low rank tensor learning. Advances in Neural Information Processing Systems, 27, 2014.
Beck, A. and Teboulle, M. fast iterative shrinkage- thresholding algorithm for linear inverse problems. SIAM Journal on Imaging Sciences, 2(1):183-202, 2009.
Bertalmio, M., Bertozzi, A. L., and Sapiro, G. Navier- stokes, fluid dynamics, and image and video inpainting. In Proceedings of the 2001 IEEE Computer Society Con- ference on Computer Vision and Pattern Recognition. CVPR 2001, volume 1, pp. I-I. IEEE, 2001.
Impact Statement
Tensor regression provides substantial advantages tradi- tional regression methods, particularly in settings involving multi-way data such as video analysis, multi-modal sig- nals, or high-dimensional biological datasets. This work aims to advance the field of tensor regression by developing more effective and scalable models tailored such com- plex data structures. While the proposed methods may have broad societal implications in various domains (e.g., health-
Bi, X., Qu, A., and Shen, X. Multilayer tensor factorization with applications to recommender systems. The Annals of Statistics, (6(66):3308-33))
2018.
Candes, E. and Recht, B. Exact matrix completion via convex optimization. Communications of the ACM, 55 (6):111-119) 2012.

---

High-Dimensionn Tensor Regression with Oracle Properties
Candes, E. and Tao, T. The Dantzig selector: Statistical estimation when p much larger than The Annals of Statistics, 35(6):231332355) 2007.
He, L.., Chen, K., Xu, W., Zhou, J., and Wang. F. Boosted sparse and low-rank tensor regression. Advances in neural information processing systems, 31, 2018.
Carroll, J. D. and Chang. J. Analysis of individual differ- ences in multidimensional scaling via an N- N-way general- ization of "eckart-young" decomposition. Psychometrika, 35(3):283-317 1970.
Hore, V., Vinucla, A., Buil, A., Knight, J., McCarthy, M. I., Small, K., and Marchini, J. Tensor decomposition for multiple-tissue gene expression experiments. Nature Ge- 2016.
netics,
Chen, H., Raskutti, G., and Yuan, M. Non-convex projected gradient descent for generalized low-rank tensor regres- sion. Journal of Machine Learning Research, 20(5):1-37, 2019.
Kolda, T. G. and Bader, B. W. Tensor decompositions and applications. SIAM Review, 51(3):455- 500, 2009.
Lee, W. S. Collaborative learning for recommender sys- tems. International Conference on Machine Leurning, volume pp. 314-321. Citeseer, 2001.
Chen, J., Scarlett, J., Ng, M.. and Liu, Z. A unified frame- work for uniform signal recovery in nonlinear generative compressed sensing. Advances in Neural Information 8252, Processing Systems, 36:8224- 26:8222-8252, 2023.
Li, D., Chen, C., Liu, W., Lu, T., Gu, N., and Chu, S. Mixture-rank matrix approximation for collaborative fil- tering. Advances in Neural Information Processing Sys- lems, 30, 2017.
Collins, M. and Cohen, S. Tensor decomposition for fast parsing with vatial-v pcfgs. Advances in Neural Information Processing Systems, 25, 2012.
Li, L. and Zhang, X. Parsimonious tensor response regres- sion. Journal of the American Statistical Association, 112 (519):113111455 2017.
Elenberg, R., Khanna, R., Dimakis, A. G., and Negahban, S. Restricted strong convexity implies weak submodular- ity. The Annals of Statistics, (6666B):353933555 2018.
Li, N. and Li, B. Tensor completion for on-board com- pression of hyperspectral images. In IEEE International Conference on Image Processing, pp. 517-520. IEEE, 2010.
Fan, J. and Li, R. Variable selection via nonconcave pe- nalized likelihood and its oracle properties. Journal of the American Statistical Association, 134456) 2001.
Li, Q.. Jiang, I.., Li, P., and Chen, H. Tensor-based learning for predicting stock movements. Proceedings of the AAAI Conference Artificial Intelligence, volume 29, 2015a.
Fan, J., Liu, H., Sun, Q. and Zhang, T. I-LAMM for sparse learning: Simultancous control of algorithmic complexity and statistical error. The Annals of Statistics, 46(2):814, 2018.
Li, X., Xu, D., Zhou, H., and Li, Tucker tensor regression and neuroimaging analysis. Statistics in Biosciences, 10 (3):520-545, 2018.
Farias, V. and Li, A. Optimal recovery of tensor slices. In Artificial Intelligence and Statistics, pp. 1394-1402. PMLR, 2017.
Li, Y.. Nan, B.. and Zhu, J. Multivariate sparse group lasso for the multivariate multiple linear regression with an arbitrary group structure. Biometrics, 71(2):354-333, 2015b.
Grippo, L. and Sciandrone, M. Nonmonotone globaliza- tion techniques for the Barzilai-Borwein gradient method. Computational Optimization and Applications. 23:143- 169, 2002.
Liu, H., Wang, L., and Zhao, T. Sparse covariance ma- trix estimation with eigenvalue constraints. Journal of Computational and Graphical Statistics, 23(2):439-459, 2014.
Gui, H., Han, J., and Gu. Q. Towards faster rates and ora- cle property for low-rank matrix estimation. In Interna- lional Conference on Machine Learning, pp. 2300-2300 PMLR. 2016.
Liu, Y., Liu, J., Long, Z., and Zhu, C. Tensor regression. Springer, 2022.
Lock, E. F. Tensor-on cnnsoor tensor regression. Journal of Com- putational and Graphical Statistics, 27(3):638-647 647, June 2018. ISSN 1537-2715.
Han, R., Willett, R., and Zhang, An optimal statisti- cal and computational framework for generalized tensor estimation. The Annals of Statisties, 50(1):1-2) 2022.
Luo, Y. and Zhang, A. R. Low-rank tensor estimation via riemannian gauss- gauss-newton: newton: Statistical optimality and second-order convergence. Journal of Machine Learning Research, 34(381):1-43 1-18, 2023.
Hao, B., Wang, B., Wang, P., Zhang, J., Yang, J., and Sun, W. W. Sparse tensor additive regression. Journal of Machine Learning Research, 22(64):1 143, 2021.
10

---

Hligh--Dienniioal Tensor Regression with Oracle Properties
Luo, Y. and Zhang, A. R. Tensor-on-tensor regression: Rie- mannian optimization, vver-parameterization,, statistical- computational gap and their interplay. The Annals of Statistics, 52(6):5583-2612, 2024.
Tucker, L. Some mathematical notes on three-mode factor analysis. Psychometrika, 31(3):279-311 1966.
Wang, D., Zheng, Y., Lian, H., and Li, G. High-dimensional vector autoregressive time series modeling via tensor de- composition. Journal of the American Statistical Associ- ation, (3337 1356, 2022.
McConnell, A. Applications of tensor analysis. Courier Corporation, 2014.
McCullagh, P. Tensor methods in statistics: Monographs on statistics and applied probability. Chapman and Hall/CRC, 2018.
Wang. D., Zheng, Y., and Li, G. High-dimensional low-rank tensor autoregressive time series modeling. Journal of Econometrics, 238(1):105544 2024.
Miao, H., Wang, A., Li, B., and Shi, J. Structural tensor- on-tensor regression with interaction effects and its ap- plication to a hot rolling process. Journal of Quality Technology, 54(5):547-566 2022.
Wang, Z., Liu, H., and Zhang, T. Optimal computational and statistical rates of convergence for sparse nonconvex learning problems. The Annals of Statistics, 42(6):2164, 2014.
Nesterov, Y. Gradient methods for minimizing composite functions. Mathematical Programming, 140(1):125-105 2013.
Wei, Q. and Zhao, Z. Large covariance matrix estimation with oracle statistical rate via intoriiatioonmmiimmizatoo IEEE Transactions on Signal Processing, 2023.
Nion, D. and Sidiropoulos, N. D. Tensor algebra and multi- dimensional harmonic retrieval signal processing for mimo radar. IEEE Transactions on Signal Processing, 58 (11):5693-57000 2010.
Xu, D. Sparse symmetric tensor regression for functional connectivity analysis. arXiv preprint vrXiv:2010.14700 2020.
Yao, Q., Kwok, J. T., Gao, F., Chen, W., and Liu, T.-Y. Effi- cient inexact proximal gradient algorithm for nonconvex problems. In International Joint Conference on Artificial Intelligence, pp. 3314, 2017.
Rabusseau, G. and Kadri, H. Low-rank regression with tensor responses. Advances in Neural Information Pro- cessing Systems, 2016.
Yu, R. and Liu, Y. Learning from multiway data: Simple and efficient tensor regression. In International Conference on Machine Learning, pp. 373-381. PMLR, 2016
Raskutti, G., Yuan, M., and Chen, H. Convex regulariza- tion for high-dimensional multiresponse tensor regression. The Annals of Statistics, 47(3)1155411844 2019.
Yu, R., Li, G., and Liu, Y. Tensor regression meets gaus- sian processes. International Conference on Artificial Intelligence and Statistics, pp. 482-490. PMLR, 2018.
Recht, B. Fazel, M., and Parrilo, P. A. Guaranteed minimum-ranil solutions of linear matrix equations via nuclear norm minimization. SIAM Review, 52(3):471- 501, 2010.
Yuan, M. and Lin, Y. Model selection and estimation in regression with grouped variables. Journal of the Royal Statistical Society Series B: Statistical Methodology, 68 (1):49-67, 67, 2006.
Russakovsky, O., Deng. J., Su, H., Krause, J., Satheesh, S., Ma, S., Huang, Z.., Karpathy, A., Khosla, A., Bernstein, M., et al. Imagenet large scale visual recognition chal- lenge. International Journal of Computer Vision, 115:
Zhang, C.-H. Nearly unbiased variable selection under minimax concave penalty. The Annals of Statistics, pp- 894 0449942,44 2010.
211-252, 2015.
Semerci, O., Hao, N., Kilmer, M. E., and Miller, E. L. Tensor-based formulation and nuclear norm regulariza- tion for multienergy computed tomography. IEEE Trans- actions Image Processing, 23(4):1678-169, 2014.
Zhang, J., Cai, Y., Wang, Z., and Wang, B. Sparse and low- rank high-order tensor regression via parallel proximal method. arXiv preprins arXiv: 12265,,, 2019.
Su, J., Byeon, W., Kossaifi, J., Huang. F., Kautz, J., and Anandkumar, A. Convolutional tensor-train Istm for spatio-temporal learning. Advances in Neural Informa- fion Processing Systems, 33:13714-13777 2020.
Zhang, T., Fu, Y., and Li, C. Hyperspectral image denoising with realistic data. In 2021 IEEE/CVF International Conference on Computer Vision, pp. 2228-2237. IEBE, 2021a.
Tibshirani, R. Regression shrinkage and selection via the lasso. Journal of the Royal Staristical Society Series B: Statistical Methodology, 58(1):267-225) 1996.
Zhang, Y., Bi, X., Tang, N., and Qu, A. Dynamic tensor recommender systems. Journal of Machine Learning Research, 22(65):1 2021b.

---

High- Dinhnsiiinnsin Tensor Regression with Oracle Properties
Zhou, H., Li, L.., and Zhu, H. Tensor regression with applica- tions in neuroimaging data analysis. Journal of the Amer- ican Statistical Association, 1085502):500-52,, 2013.
12

---

High- High--iimensionl Tensor Regression with Oracle Properties
Supplementary Materials for "Hiigh-Dimensional Tensor Regression with Oracle Properties"
Content of Appendix
A More Low-Dimensional Structures
14
A.1 Fiber-wise sparsity
14
A.2 Slice-wise sparsity
14
A.3 Slice-wise low-rankness
15
Proofs the Theoretical Results
15
B.1 Proof of Theorem
15
B.2 Proof of Theorem 7
22
B.3 Proof of Theorem 8
23
B.4 Proof of Theorem
24
B.5 Proof of Theorem
35
Complementary Experimental Results
36
C.1 Synthetic Data
36
C.2 Real-world Data
40
13

---

Tligh-Dimensional Tensor Regression with Oracle Properties
A. More Low-Dimensional Structures
In many applications, tensor coefficients exhibit structures beyond the aforementioned. For instance, an entire fiber or slice of tensor might be zero, or a slice might be low-rank (Li et al., 2015b; Raskutti et al., 2019; Chen al., 2019). Below, we introduce several penalties that promote these more intricate structures, followed by a corresponding theoretical analysis.
We now provide oracle-style guarantees for estimators employing these penalties. Our analysis parallels that of Section 4 and leverages similar assumptions.
A.1. Fiber-wise sparsity
Consider unfolding the tensor A along mode-(k), the fiber-wise spansity regularizer is defined as
Σ pA (I[Av)).lb)
RS (A)
where the [X].; . denotes the j-th column of X and |I-I|2 is the vector 2-norm. This penalty encourages entire fibers to be
zero analogous the group Lasso (Yuan Lin, 2006)
The oracle rate refers the statistical convergence rate of the fiber-wise sparse oracle estimator, which is assumed know
the true support set S3 in advance, where S3 = i [A0)
. |1. + 0). Specifically, the fiber-wise sparse oracle estimator
is defined as
A, arg
min A:A
C ((A).
V d and [A satisfies the Theorem (Fiber-wise sparsity). Suppose that Assumptions 1 1-4 hold. 1f > A= X= condition that
min IESs
z0,
the estimator A to problem (2) satisfies
SSa Ss de
A
A
2 F
A.2. Slice-wise sparsity
In other scenarios, one may expect entire tensor slices to be zero (Raskutti et al., 2019). To promote such slice-wise sparsity, we introduce the slice-wise sparsity regularizer, defined as
RA (A) = negg
Σ pa (((Al))
The oracle rate refers to the statistical convergence rate of the slice- wise sparse oracle estimator, which assumed to know the support set Sa in advance, where S4 {1 [A and 1|: * +0}. Specifically, the slice-wise sparse oracle estimator is defined as
Ac arg min A.As,
L(A).
Theorem 8 (Slice-wise sparsity). Suppose that Assumptions 1-4 hold. 1f u > S, X V condition that
d,d2 and [Acsko] .nd satisfies the
min LES,
|A
24,
the estimator A to problem (2) satisfies
À A'l
Sa d
I F
14

---

High- Dinennionn Tensor Regression with Oracle Properties
A.3. Slice-wise low-rankness
Beyond mode-wise low-rankness, some applications benefit from low-rankness within individual slices (Lock, 2018; Raskutti al., 2019). To capture and exploit this structure, introduce the slice- wise low-rankness regularizer, defined
= Ummaa dm min(d,de. -eess de
PA (5. (1A6))))...))
RA (A)
The oracle rate refers to the statistical convergence rate of the slice-wise low-rank oracle estimator, which is assumed to know the true rank S5 $ estimator is defined as
as
[Tx
((Ais)))..
)
* 0 advance. Specifically, the slice-wise low-rank oracle
A arg min C(A). A.As,
Theorem
9
(Siece-wise low-rankness).
Suppose that Assumptions 1~4 hold.
If P Ç x
x (E(L))
([AciA)..) satisfies the condition
S
and os
(5.K) ..d sddlep
a. (145sN)),. m11) )
2 2
[[x [xx(((())) (E(Am))
}
2
nppe
(1.b) 1..
the estimator to problem (2) satisfies
T(k / Ss
J A'
A
<
n
where T(.A) = max [1, (IE* (E))
B. Proofs of the Theoretical Results
B.1. Proof Theorem
We begin by demonstrating that the entry-wise sparsity regularizer can be reformulated as the sum of the & penalty and a concave part. Specifically, we have
RA (A) =
Px (Aig....ix) = X ||at,.... 11, + Qx (A),
1==1
ty=1
where (dlp....yll = -qp-1 Edy
lai is the E, norm and Qx (A) = Em Eay aa (ai
Lemma 1. Under Assumptions 2 and 3, the loss function E (AA satisfies the restricted strong convexity
[AA) 2 ((A) + VĪ A)
A
All2
and the restricted smoothness
Σ(A') Sī(A) (Vū(A) (VŽ())))(((())LL.)))))))) 4' A)+ + |
L All2
Proof. Recall that Q, (4) represents the concave component of the non-convex penalty R, (a:,....ix) imply- ing that (A) is convex. Specifically, 2x (A) can be expressed sum over its entries 2x (A) E Eiy= aX (at,...ww) where au (a.,...,,x) satisfies the third regularity condition specified in Assumption 1. From this assumption, we have
-C (dlg...ii Oig...i )² S (o' (d/ 4 (ai,...ii)) (a/g.mxx aijp..ix) 0.
15

---

High- Bighh-iimensionl Tensor Regression with Oracle Properties
By aggregating over all entries, deduce that the convex function -2x (A) satisfies the following inequality
(VV -
.A' A S A - All
(6)
- ( -
A'' A 2 0.
Inequalities (6) and (7) correspond to the definitions of RSC and RSM for the function -2x (A), respectively. Specifically, they imply that -2 (A) is both (-smooth and 0-strongly convex. Consequently, we have the following
- Ox (x) - -2x (VQ) (4) A'' - A) + 2 || - All
- Qx (4) 2 2x (A) - (VQ) (4) A'' A)
For the loss function C (.A), applying Taylor's theorem and the mean value theorem yields
L (AA)) = (A) + (VL A'' - A) < + (1 -1- A) (A' A) @ (A' - -A))
for some / E [0, 1]. Here, @ denotes the Kronecker product. Given two tensors A, A' E Rdix-xxx their Kronecker
product results a tensor A" of dimension (did1) x x... (dydv). Each entry di iafa..aaa.... is defined as ai da X ji.
Under Assumptions 2 and 3, we have
(IB)-(AA) (VL(A), B- - A)
+ Ale
2
T
( C S (VL B A) 1 AlE:
Recall that (AA) (AA) + 2x (A). Thus, we obtain
Σ(A) 2Ī(A) + (VĪ(A), A' -A)
--
| - Alp
2
L
[AA) SÉ(A) + (VĪ(A) A' -A) A AA
Lemma Suppose there exists an integer $1 > C\Si), where is a constant, and that A satisfies || S 81. (A) < 3 and IIVC (la S X/8, where || Ilmax denotes the maximal element of the tensor. Under Assumptions 2
and 3, A satisfies
21/8
I A 1r S μ X
Proof. Given that || A
lb S ड. and
A
0. it follows that [I(A
)5,11 S 31. Based on Lemma 1, we can
3,
derive the following inequalities
(VÊ(A), A)
/ リー6 | A 1
((A)2Σ(A)) + (VŽ(A)),A Ar)
/ 4ー6 I 1
Adding (8) and (9), we obtain
(VΣ(A),A'AA) 2) A- A')+ + (4-)) A
(10)
16

---

High- Bighh-iimeniionll Tensor Regression with Oracle Properties
Let E 0 | denote the sub-gradient and S be set. According to the Karush-Kuhn-Tuccke (KKT) condition, we have
(ā( + GG =0
Determining the optimal solution is challenging, therefore, we introduce a measure of sub-optimality
AA-A (4 (A-A', (A=A.VE(4)+
1
w (A) min max A'ES
9 e'E
We define our algorithm's 's stopping criterion as w (A) S s. Consequently, the sub-optimality can be expressed as
1
(A (4-A, +
(A)
max E'
A A'1,, ||,
Adding A (A- g") to the both sides of (10), obtain
A - A" ( > A A'. (A-A).VE(A) +
V +m
Al AI(AA-A''G)
Since A' E ES, we have
[A-Ali (A- (A-A A', v2(A)) +X0) ) (A()
A - A'll, A", VĪ + +10 S max
A'ES
Recall that we assume (A) S 2/. we obtain
A (A- A4 +1) s0(A)
11 A:ll
(11)
2
Combining (10) and (11), we obtain
X 2 -1- 2 A A', A', V((A)) + (u ( -) + A'
For term 1, separating the support of A - A" into S and Si, we have
A--;,V2(A"))
A A vā(A4) + (A A (2,)
2 |4 A (AA) mlx (A- A", ve, (Ar))
= [(4 13s, |1, |/V.C (A4)) - 11(A - A')s, Il, | (Ar)
+ (44 - AT)s, (VQ, (s) + ((A A')s, (VQx (A)))s)
2 (A A') ||VC(4 Il Imaxx
I|AA AA'ss, 11 IVC (A4")
(A A s1, /| (A4*) I Iaax
For term II, separating the support of A - A into S₁ and 3,, we have
A A' 9) => ((A A' A')s, ,s))), + } (AA--AAA - A')
2 A - A')s, |1|
|9 Imax + > +A(A5, 95,)
2 --AI(AA- s7),, +X Σ
lais..it
(H.....NNES.)
=-AI(AA - A'), 11, + +A|
- -A')²ll
17

---

High-Diin Dishnsinnnssin Tensor Regression with Oracle Properties
Thus, we obtain
A 2
A
A'll1 2 - II(A A') ST
[V mIm || A:)s 11, IVC (A*) max
(12)
- -I[AA- (IIA-A') s'),, || | (44 Il + (4-C) () || - All
A/I(A- A'),, |l + - -4),,1111--
We separate the left-hand side of (12) as
マタ2 |-- A'll = - II - A')s. 11 + 2 I[(A- - A')
A
A
A
2
Rearranging the terms, obtain
(u 1)-5) - All + (2 (24(4"))max) ) II(A A')
31 2
<
+ IIVC mmax + +I/V2x (44) 11 Imaxx |I( A1),111
Recall that VLL (A'))
maxx S
8. we have
A
3\
(4-5)
All s
+
[|VC(A')]
max + 1V2 (A4)l
II(A (A4-A))s. 11
2
3X
X
2) |I(A AA')s,
+
2
8
21
V (A - A')s.Il
8
21X
VIS
|A AAI
Given that a - <> O, we have
21/8
|| A S 1
X S
Lemma 3. Consider the regularization parameter À and assume that the derivative of the non-convex penalty satisfies
p'a (ais....i)) = O whenever (N 2 V for some >> O. Let St U U = Sy. For indices (f1...., in) E S C S₁, we assume o.... 2 v, for indices (i.,... iN) Sit C S₁, we assume [ S v. Under Assumptions 2-4, we derive the following bound
A
Arl
1
| + 11 "u-c" --
3
> SP-
Proof. Define the sub-gradients 9* @ and & E a ||
Note that A satisfies the optimality condition that w (A)
s 0, we have
max A'ES {(2
A', Vī (A)
+ +10)} so.
(13)
18

---

High- Bigeh-Dimensioal Tensor Regression with Oracle Properties
Given that As
S 81. since (a-
S S1, according to Lemma 1, we obtain
E (2) 2[(A2) (VÊ .
1 A
(14)
E(A) (A)
(Vē (A) A - -A)
||A Ali
By the convexity & norm, we have
2 ||j|| > + -1(AAA1'"))
A | S -A(A--A.)))
(15)
Adding (14) ~ (15), we obtaín
VL (L() + VQx (a 0g*, A A1) +
-VV A
A' +)+
+ (1 (w- (m- A (16)
(i)
From the optimality condition (13), we have
(Vī + A A'
A) S
+ +10) 00,
max A'ES
which implies the term (1i) in (16) is non-negative. Consequently, arrange (16) obtain
( - <) - |
VC(A + vys (A)) + + , AA-A))
IN (VC
(4- 2) (e..d
min Ar1l
Σ
VQx (A)+ + + (
(17)
1y=1
We proceed by decomposing the summation on the right-hand side of (17) into three distinct parts
(i,,..., iN) Si.
(i1,...,4n) E S..
(i,....iin) E S!I
Here, Si = ((i..... iN) |as,..ix 2 v), St ((1....,, ix) | < v), and v > 0 is defined in Assumption
(i) For any index (i.,..., iN)) Si, the regularity condition yields
VQaaj =44 (at,..x) d (0), for j E ST.
Assuming that [V'' [Iax) < À g it follows that
d S X.
max (N).
(VL(A))) (
I VVLC
Therefore,
max (3)
(VL + Qx (A4)). ) 5X
19

---

High- Diahh-Dimeniinal Tensor Regression with Oracle Properties
Moreover, since g" 0 || it holds that Xg IN....A ( Consequently, for each (6.....,nn) E S1, select g i such that
(VV + VQx (A4)),...s -ggg... 0.
This implies
min ArlliA
([V VQ (4*) +2')))+++))"..... (..) =0, for (b (fI...i)) E S1.
Therefore, we obtain
Σ (VL + (4') Xg") x (A' -A) (I..EN 1o. 0.
min gred A
(18)
(ii) For indices (d1....,dN) iN) E S we have | 2 v. Given that R(A) = 1|/ + Q2 (Di our assumption
on R (.A) ensures that
This leads
min A-E
Σ
(VL(A)) VOA
Vo (A*) + ) ( (5-AA -A) $1...
(A'-
(H...N)EES)
Σ
IVVL (Ar))) (Am)))
(4:-2) i...iN
(VL (A
I|: |
A - Al
(iii) For indices (iN)...,i)) iN) SHi we have |07 (a] Given that VL (A')
max Ilmxxx S 8.
we have
(VL(A) )) S (VV((A))) )),,.. max S $1/8.
(1 (ESH
Meanwhile, we have
(VQx (A'))
max (i...iNJES
10 (a))..x) S max 14's (o ) S 2,
(t...in)ES))
Additionally, since 9 0 || it follows that I S Therefore, for each (d1,..., iN) E E!! obtain
(VVC(A)++VN + VQ: xg*) 1ynnn
max (Hi))))ES))
(VC(A)), sqnnee max ((
(VQx (A4)), s...dy
< 3A,
which implies
Σ
(VL( + VQ. (A") + >g* >g"))...... )nn..wl (A' A)
min G'E ||
ig...N
(A' - =)
31
AAt A
3\
S |
(4' A
<3\ |
<33 S
A A
(19)
20

---

High- Diahnsinnnn Tensor Regression with Oracle Properties
Substituting the bounds from (18) to (19) into the right-hand side of (17), obtain
A' A|p S --s (I(VE(4)))))) + +3 +3 |/
Lemma For least-squares regression with sub-Gaussian noise, we assume that the columns of X are normalized
in such a way that maxje (l....( Xd2 xxxxx & 1/
I vn, where X (vec (x(1)) p.... vec (xim)) If
log do- dy) dv) then we have
S n
(A 2
Proof. We begin by establishing an upper bound on the probability that the maximum entry of the gradient (I)V 2 } where (A) = x E and & (vec (E()) vec (g(n))
For dogid using the union bound, we obtain
)
2)
c colog d/n
VLC max
SP
8
72
8
di x xdxxx..xdx
Σ
P
(x)
na
E)
2
Cv
log d/n
8
Let's define 0%
= (x. where & is composite coordinate. Since E, is sub-Gaussian (0. (0. it follows that for any
to >
1 1//2/22
2
(exp {Lo02} + exp {-100x3) = 2exp
X
n
Taking to
2n
yields that
| 202
nt2
P (1ON 2 1) S 2exp
2 ||1x(x1
x
2
172
Further taking results
1 8
xxxxxxxxx((1x)) -22((1288²)
IP VL (A) (A) max
< 2 (di ... ....xdyy xd)
Applying the Hanson-Wright inequality yields that
P ([
(A, E))
EE, BA
E(E, (A.E))
(E(E)
}
exp
Cmmnn
( |All
172
A
where C' a universal constant.
21

---

then we High-Diimeniooal Tensor Regression with Oracle Properties
Combining the above inequalities, we have
2E ε, X
S
VVL (F
I v
v2
(20)
n
Building upon Lemma through Lemma 4, we derive Theorem 5.
B.2. Proof of Theorem 7
We begin by demonstrating that the fiber-wise sparsity regularizer can be reformulated as the sum of the & penalty and concave part. Specifically, we have:
1l d2
Im dj Σ
Rx (A) Σ PA (1|AAr)),.11))==
> 1|A4 + + (1(AA))).lll))
where QA ([AAm))))l) Elliend aa (1((A()))ell)
Lemma 5. Under Assumptions 2 and 3, the loss function E (AA satisfies the restricted strong convexity
[AA) 2 ((A) + (Vī(A), ( A)
1 C |
All
and the restricted smoothness
SQ)))(()) (VC(A) A' A) + - || All
Proof. Since the proof closely mirrors that of Lemma 1, it is omitted here for brevity.
Lemma 6. Suppose there exists integer 83 > C1S)), where C is constant, and thut A satisfies ||4
83
min (A) S 2 where w (A) geal [wwl.
V E (A) Xg
and V (I S 2/8. Under Assumptions max
and 3. A satisfies
21/8
||A A S J
X Sz.
Proof. We omit the proof here for brevity, as it closely mirrors that of Lemma 2.
Lemma 7. Consider the regularization parameter À and assume that the derivative of the non-convex penalty satisfies
p' = whenever
A | 2 v some some v > 0. Let Shush! = S3. For indices (i1,...,iNN) E S C S3,
assume min Aix) 2 u, and for indices .i.....,) S. C Sz. assume min 2,4, derive the following bound:
(Ai)) S Under Assumptions
A [ A-
S - (VC(4))) + # -
1
3 A 1
Proof. For brevity, we omit the proof here, as it closely resembles that Lemma 3.
Lemma 8. For least- least-squares regression with nab-Gaussian noise, we assume that the columns of X are normalized such
X S nn, where X = (voc (x()) vec (xx(1)) If X x v logga du
way that have
S
IVC (AA")l S
n
Proof. For brevity, the proof omitted here closely follows the methodology established Lemma
22

---

Bighh-Dimensionll Tensor Regression with Oracle Properties
B.3. Proof of Theorem
We begin by demonstrating that the fiber-wise sparsity regularizer can be reformulated as the sum of the & penalty and concave part. Specifically, we have:
I 1 d,
IT Nep d, Σ
RA (A)
Σ
Pa (11AA
A + + (1(AA)s))millll))
i=1
==1
where 2x (AAAAAA) .lee Monsi as (I(AAOM))....ll)
Lemma 9. Under Assumptions 2 and 3, the loss function E (AA) satisfies the restricted strong convexity
2 + (VE()) A) + 1 P-5
|
All
2
and the restricted smoothness
{(A) SC(A) (VÍ(A), (VŽ(A)) A)+
1 || All
Proof. The proof can be demonstrated similarly to the proof in Lemma 1. Hence, we omit it here.
Lemma 10. Suppose there exists an integer 34 > C. (1)) where C is constant, and that A satisfies || 1s S 34 w (A) < Σ.. where w (A) =
[1 (A) + Ag } and |V (AA)) Ilmxx S X/8. Under Assumptions 2
min (eal[A(s)))
and 3, A satisfies
21/8 Sa S4
IAA A
Proof. For the sake of brevity, omit the proof here, it closely follows that of Lemma 2.
Lemma 11. Consider the regularization parameter X and assume that the derivative of the non-convex penalty satisfies
P (([(()))))))))))) =
whenever
Agjii.... 2 v for some v > O. Let St U S = S4. For indices (ii,..., iN)
2 v, for indices (i....., iN) (x)) E SH C Sg, we assume min v. Under Assumptions 2-4, we derive the following bound:
S C S4, we assume min
[AG -
[A ... S
A - A-1
1 S μ ((V((44)))
3 μ A SHI-
+
Proof. For brevity, we omit the proof here, as it closely resembles that of Lemma 3.
Lemma 12. For scasres regression with sub- Gaussian noise, we assume that the columns of X are normalized in
such way that max X ||2 VT, where ă = (vec (x(1)) ..... vec (x(1i)) If X log(d,ds)
then we have
Sal
VVC (AA))l S
Proof. For brevity, the proof omitted here closely follows the methodology established Lemma
23

---

High- Dimensionn Tensor Regression with Oracle Properties
B.4. Proof of Theorem 6
The proposed mode-wise low-rankness penalty can be reformulated the sum of scaled norm and concave function. Specifically, we have
min minn(ds.I1,,t 14 d
RA (A)
Σ
pa (A(e)) = Ax) + Qx (A()).
where (A()) denotes the i-th singular value of the mode-(k) unfolding A(a)- For the estimation problem, we define
Σ(A) L 2x (Ax)),
where 2x (A())) (A())) -. (n.Tl ds} aA (o, (A())).
Based the restrict strongly convexity C (-) Assumption and the parameter for regularity condition in Assumption 1, 4 > 5, we have the restrict strongly convexity (.).
Besides, for the RSC and RSM assumption, we define the following cone of directions
C { IR dR ||| (B))lmu I [[y (B) nuc)
Lemma 13. Under Assumption 2, if B C.we have
Ë ()+BB) B) (V
1-5
B |
2
Proof. Based on Assumption 2, we have
L (((A)
|BF
(21)
Moreover, considering the singular values of the unfolded matrices A(k) and (1)). we obtain
-SS & (a: (A())) - 9 (a: [AA ( B3 11()))
o, (A())) - oi (IA+ + Bl(x)
which is similar to the proof for Lemma 1. This inequality leads to
- (-VO, (1++ B(n) s B
This inequality characterizes the smoothness of -Q(-), which is equivalent to
2x (IA+ (Ba) Qx (A()) (a)) 2 2x (A() + (vQx (A(a) B(x)
C |
(22)
Noting that the Frobenius norm satisfies |B | = BB Let A' = A + B, adding (21) and (22), we have
= 2 (A') + Qx (A(t))
2 (AA) + 2 (A(A) + (VC BB + (V2x (A(s)) BBb)) (A()))B())) + 1 C
|B8
2
=Ā(A) + (AA) B) + (VQ) (Ax)) B(x)) + 11 -
| 1
|
2
24

---

High- Dineesiimnniin Tensor Regression with Oracle Properties
(e)( 2n
Lemma 14. Under Assumption 2, if p. > { (and the regularization parameter À
we have
| (Am)
55 my (Am) A(s)) )
Proof. By Lemma 13, we have
ī -())
(ve (Ais) | AA'
(23)
(k)
We proceed to bound the right-hand side of inequality (23). By decomposing the inner products using the projections onto two orthogonal subspaces, we have
VL(A') A 4:) (vex (A(x)) [A A.) (10)
[VC (A (M)) +v2 (Ais). _1I= 115 (1a-A*) 111m)) (VVC(A)))ss +Vex (Ais) F1ا (12- A.) (1)
x1((m))
2 --|0r (VL(A)))() Uee (4ix)) [I, (la-A))
(A-A1 _A.] (N)) muc
(24)
- |ر (IVC(A)))()) + vex (Aie))
1z (12- AA) (k)) 10n))
(25)
For (24), due to A 2 21 x (E)) sp
we see that
[VV
S/22. According to Assumption we have
|1 (VC(()))()) ver (4ix))
3
S
(26)
2
For (25), since 1p (Ais)) = O, we obtain
درا (IVC(A)))()) vex (Aix))
1
(27)
2x
Combine (26) and (27), we have
(VC(Ar), A- - A+) +++/oe (V2x (Ais)) [ [A-A*) (10)
- マタ
3
||I i u5
A-) (=(2))
1 2-1
(1.2 A-' (12))
براا
Moreover, noting that X bining with (23), we obtain
A
-
A || ||
2 2
2-1/us= [1> (1a- A-)
11)))
)1
+ 11ا1 ( ([a-A') ( (10)
|1
and com-
(3))
VL(A') + +Vx (Ams)
AAA + |
À -
-A11A" Ilane
;2 |1 II (1ā- A")
(12- A-)
(16) 11))...
x 2
(28)
muc
Since A is the global minimizer of the general estimator (2) and given that μ Ç, it follows that
A ( + 2|| nuc
-(AA) A ||A Ianc < 0.
(29)
25

---

High- Digh--iimeniooal Tensor Regression with Oracle Properties
Substituting (23) and (29) into (28), we obtain
2 در ([2-Ar)
A+ (8)) nuc
)1
11F (12- - AA (10) \ S 2
Since À > 0, we obtain
|مر (1 A 10a-)r10)).1 nuc
5 |1=
(1a-
Lemma 15. Considering the mode-wise low-rankness regularizer, under Assumptions 1~1,for 1~11 the estimated parameter tensor A and the true parameter tensor A', we have
2 A
1 S p p-S
[VISH]
|1= (1Z((4))))),,, + +3V S
where St and S are subsets of the support set of St. The set S include all indices i E S which satisfy o; (Ais) 2 v.
and S includes all indices with (Aiss) < v.
Proof. Since || llaue is convex, we have
A
((s)
A A1s)
+1 +2(12-
(12- A+
G*
(30)
mac
mur
(N)
>
A(x)
2 X A (1))
+ (( 21 (8)
G
(31)
muc
nuc
where G* (x)lAx and G E 0|| (A)) Iluuc- From (30) and (31), we have
A I + +>|1ATo))
2 x ATo
+ +
+ +x (1 A+
(0)) G" +1 (14r-
A
GG)
muc
(K)
(b)
This equals to
0 z x Ar) (0))
: G"
+1 (14
2
(32)
(1)
Moreover, according to Lemma 13, we have
7 (())2ū(A)) E(A') + (VC(A)) ( (--)4) (
vex (Ais): [Ā [Ā-A*) A- (k)
/ M- C
A A
(33)
Z(AT) 2 2 (²)) + (VC (?) A 2)+ +//e, (A(1)) [A' [A--AA] -À
-2)+
+ / リー C
A"
21
(34)
(8)
2
F
Summing (32), (33), and (34), we have
02 (L(A') --A + (VC (A) .A' 22) + +(u- 13-AA-111
4+)
A..
+ ((Tee (Ais) + +X",, [ [Ā- A- (k)
)
+ (ve, (A1t)) +XG, Air) AAi) -A1))
26

---

High- Bigh--imensional Tensor Regression with Oracle Properties
Since A is the solution to the estimation problem and A satisfies the optimality condition, for any A' E Rd xedxx, it holds that
A 1²) (ve, (Awx) + +Λद; [A- {(vC(()
- A (8) )
max A'
10,
which implies
VL (A) AA -4-24 -2) + (ve, (A(s) + A2 [A - A) (k) 20.
.A'-
Since (VC(A"), A 2) - ( VV (00) [ Al
we have
(2)
(n (4- 11 1--00 A < ( (V M00) |A-- -2- (10) + Vex (Ais) + +AG", [A--
-A - )
< 11 [|CC(4)]()) + V (dit) + +XG*
114--2] -À (m)
II (VE((A))()) + +Ve (Ais) + +AG*)
[A- Al
(35)
(w))
We have defined o; AA))) as the i-th singular value of matrix Alay. With regard to the magnitudes of the singular values of A (k). we can decompose (35) into three parts:
E S that o; (Aix)) 2 v,
E Shi that 2 0z (Ais) > 0,
i E St that oi (Ait)) 0.
(i) For i S that oi (Atx)) 2 v. define a subspace of F associated with S as follows
Fst (U" V*) = [W] row (W) C Vi. (W) c U...
where Vi and Uf is the matrix with the i-th row of Vi and Uf with i E Sl.
Recall that R^ (Ais) = A7x) buc
+ +2x (Aix) we have
VRA (Atis)
v2 (Aiss) > (UiViT + +zi)
where - =--'1=
)())) Since |zi| < 1 and Zi E Fs₁ which satisfies the condition of W* to be
sub-gradient of ||| | Projecting RA (Ais) into the subspace Fsy. we have
TTss (VRA (4ix) ) =Tr (V2, (Ais) + AU + +Zi))
U14 (Ef) V + +AUGViT
UT (ET) V
where J is an identity matrix with the size min (d., II [dd } and (0's (Ei) + AT)) is diagonal matrix that for i $ S1. the i-th entry on the diagonal equals 0, i.c. 19 (27)) + A11k O. and for all E Sl. we have
1 (Ei) + A1I 4 (a. (Ais) +x= p' (a. (Ais) )) 0.
27

---

High- Mighh-iimeniooal Tensor Regression with Oracle Properties
The last equality is derived from fact that E S satisfies Assumption 1. 1, (1) -0. Therefore, we have d (Ei) + AI1 0.
which indicates that IT (VRA
(Ais)) O. For G =U'Vi + +Zi E a A(x²) we have
muc.
F/ [VC(A))) + +GG" + V2, (dix) ) A Al (K)
[1758 [VL + +VRA (dis)] [ -2] (R)
[At
(IITA) ((T((4))()) .Trss (14-- -A] (K)) 1(m))
(K))
- | (IVC
1|. up
||اa| (L - A (A)) Inn
(1)
where the last inequality is derived from the Hölder inequality. For [ias (14-- A (k)
)
from the properties of
muc
projection on to the subspace Fs!, we have
Tr3t (14- -4--2 2] (1))
)
Il.. - VISS
-iss (14-- - Al )11))
11
F
[ - Al 111 = =Vish/
A -/111
(k)
We obtain the second inequality from that the rank of the matrix MF (14--2)
S $|S
Thus, we have
(8))
11 [VC(A + +Vxx (Aix) + +XG)|
[-A'- A1
(K)
s s/S |hl |ras (I4- A1 (6)
) ||-- Al F
(36)
(ii) For E S7, v 2 oi (Aix)
> 0, define a subspace of F associated with S as follows
Fsy (U* V*)= = { row (W) C Vin col (W) C Ui))
where Vi and Ui the matrix with the i-th row of U and V with i E S!' Obviously, for all W, the following decomposition holds
TT (W) Tr (W) + Fy sn (W).
In addition, since U*, V are unitary matrices, for subspace Fst and Fsg we have the complementary subspace F FL thus we have
Fsz C FFyy. and Fsy C Fsy.
Similar to analysis in (i) on St. we have
Hf (VRx (Ait)) Unihk (Eii) Vi
where of (En) is a diagonal matrix that [4A (Eii) = 0 for i & sHl and for all i E s/,
[os (a. Aiis) )){, S.
10 [ (Ei)la
28

---

High- Dimensioon Tensor Regression with Oracle Properties
Since oi (Ais) S v, and qu (.) satisties the regularity Assumption 1, lg (t) - X. Therefore
Ial (Vex (Ais)))l,
= max M44 (Ei S A. iES!
Mcanwhile, because of the fact that Fsp c Fs. we have
|iraa (AG') - Mr (AUttiT)]!.)
(37)
|
Since
=
1, we have
sp
|> (JU
| =
(38)
up
Thus, from (37) and (38), we have
1y (AG*)
SA
(39)
sp
Additionally, due to the fact that | (IVE (A))(()))))
)1,p 5 <|IV 113) 11pp
X, which indicates that
1y (IVL(A))(() + vex (Ais) + +XG*) :..").((' ((4-- - 1000
[Iaa (VZC(A))(2) |A- (K) ) + [IFay (Ve, (Aizs) [A-- - 2 (k)
+(UEFa)
)am) + +(ura By (AG"), [A-- A] (M))
(Juss (VC((A))s)) (8) )1 + +| (Vex (Ais) )11b + || (GG3)l )
||1ur ((4--2) (1)
1 )
nuc
where the last inequality is derived from the Hölder inequality. Since we have obtained the bound for each term, as in (38)
and (39), we have
To (IVL(A))(s)) + vQx (Aix) + +XG+) [A' [A--2] -a
([4: -A
3A
IT s4
(8)
(A)
S 333 S
A - Al.
(A6)
= ST
A - A
(40)
where the second inequality utilizes the fact that rank Terae (AA- -
) )) - |
(1))
(iii) For i E Si, which correspond to the projector H since of (Ilzu (Ais)) = 0.
Based on Assumption 1, 9사 (0) = 2 (0) = 0. We have that VQl (Ais)))) = UU' (E:) V where Σ. E R is diagonal matrix and r = min(dR, [Ij ). Now we have
1 [exx (Atx)) (I. -.U. U2ah (E:) V. I - -vv.
-U: -5--U -U:) A (E:) (veT_ V
0,
29

---

High- Bighh-Dimensionl Tensor Regression with Oracle Properties
where I is the identity matrix. Meanwhile, since
x (x) llgp up
| (VL(A))) 11 5 <|
S A.
n
For -->-11111- (VL(A')) and G" = U.V. + z E 20|A nnce
and
we have
1| [VL(A + +XG* =1 ((TC(4))( + +AZe =0,
which implies that
Tyy (IVL(4))() + + + vQi (Aix)): [4' - -7] (N)
= -(00 0.0 A7 (k)
)
=0.
(41)
Adding (36), (40) and (41), which indicate that
9) A 1-AA- 1-
(p
11= ( (A ((()) + VOX (Ais) +AG*) A4' - A
-A'-AA
11---(1 -211 +32 / (A)
$ S |Ts (VL )l) 1|----
A A
= A A
|IS !وك (VC 11 +31 |/ls
Thus, we have
A A-
A
||
- μ p-c
1
VISS |eا (VZ((A))() 1l, + +32y S
Lemma 16. Suppose A E IRdxx-xxxx with rank of each mode-(k) unfolding |Sa]. Then the error bound between the oracle estimator A and the true A satisfies
2 S1l
(IVC(A'))
A A01
- 2 -A
(42)
(60
4
Proof. Let B' 1= 20 - A". According to the general estimator (2) and the definition of the adjoint operator X(-), we can express the difference in loss as follows
1 <(²0) L(A))
1
y0())
xx0) (A* + +B))
2
1
(0(1) x(() (A
2n
2n
1
E(()) x(() (B')
1
(1)
2n
2n.
1=1
1
x
(B')
1
xx (8) BB"
2n
sp
Given that A minimizes 2(-) over the subspace F and Aln) E F, we have
c - - -C S..
30

---

High- Dishhsiinnniin Tensor Regression with Oracle Properties
Thus, it follows that
| (B (B)) 2n
(1)) S (x (E) (B).
(43)
By the RSC condition 2, we know that
- (VC(A), (VE(A) + |B3
Applying this to B',
1
BB (B') L(A") - (VL(A'),
1
(1)) 1,
1 n
x (B3) )1
(x" BB) - (V)())))
(44)
2n
Substituting (43) into (44) gives
B
2
1
x
(B'))
2 (19)
1
(2)
2n
ap
n
Therefore, we have
2(1 [I (1))m))) .B) .B8
||x (Ix* (E))())
-||8')
B' 5
nps
ni
Using the fact that rank (B') = =|Sa/. we have
BB(x) Imme S <vis
B'x)|p
Thus, it follows that
2 S4
[T (x (E) (E)
B' 1
B(()) 1| s
n/
Recalling that VC(A') = (E)
we conclude
2 Sa
II [3E* (81())
2//S.a| (()))))))))))))))
B'(o)
ripu
1
Thus, since B = B3 F. we have the desired error bound
2.
St
11 (VL((A))()))
A
A'
=
|BB
B
5
F
Then, we prove Theorem
31

---

High- Bigeh-Dimensonal Tensor Regression with Oracle Properties
Proof. Suppose G e (A(s) )l
since A satisfies the optimality condition, for any A' E IR xxx it holds that
Iuc"
(VC (²) A 11)
+ v2x (Aus) + +AG, [ - A')
)
max
- 0.
(45)
(k)
GO E0 c0/
In the following, we will show some
satisfy that
muc
{{[TC (2 (M)
+ vQx (Ait) +xO", [ - A 00.
(R) )
max A'
(46)
Recall that (AA) ((A) + 2x (A(s)). Projecting the components of the inner product of the LHS in (46) into two complementary spaces F and F=, we have the following decomposition
VC (2 (K)
+ vQxx (Am) + +xG [29 A (
(1))
VC (7))
+ vex (A(a) + +XQO .1l5 (120 -A')
(K)
(16)
VC (49) + Vel (A()) + GG .I
,Is (120- AA1
(47)
(k)
For Term P.. By applying Weyl's inequality for singular values, we obtain
max o. (Ais)) - oi (Ais) - Aiss - 1.
Further, from the properties of the Frobenius norm, we have
| A 2 = 1A:22
(k)
From Lemma 16, the estimation error AA' A 0
yields
2 Sal
|x (E))(s)
max 05 (Ais)) oi (A1))
np
where [S,] denotes the rank of the unfolded matrix Ana-- Utilizing the weak condition of the singular values, we find
min
(Aix)) v
2/18a np
[ (E))(s)
Applying the triangle inequality, we derive
(A(s))
min of AA()) 5: (Ais) oi (Aïx))
min xe5,
iESs
max ES
a;
(A(1)) as (Aiss)
+ min
0i
(A)))
ies
2 2//Sal
||:0 (E) (E)](a)
+ +シ+
2 Sa
| (8)
(k)
n/
32

---

High- Bighh-iimensionl Tensor Regression with Oracle Properties
Considering the definition of oracle estimator A E F, which implies the tensor rank of each mode-(k) unfold-
O ing rank (k)) = |Sa]. And we have the singular value decomposition A(a) = U V*T. Since R2 (A(s))
||| ||) 2x (A()), for 20 Ft,
1, and Es. Sal is diagonal matrix, where
11f Φ E ) & (SS.) Based on the definition of VQ (+) and Ø I1-Ilmoo we have
TT (VR (Am)) = TT (ex (Aia) + 10
| A
yl UU'ys Σ V +AU + +12²)
= U* (
[S. + tAls.) V
(48)
where the second equality in (48) is to simply project each component into the subspace F. Is, is the identity matrix and S. E [[Su|x|Sa|, Since pa (1) = ςx (1) + A [4], we have y (1) = 9 (1) + At for all I. >> 0. Consider the diagonal matrix
dx (s.) + s.s. we have the i-th (i E S₃) entry on the diagonal that
h ES.) + Nss. = 1/ (0: (A(a)) + +A= +>=p's (os (A()))
Since PA ()) satisfies the regularity condition (ii)) in Assumption that pb, (t) = 0 for all 2 vs we have p/s (or (A(i)) =
0
for i E S4, due to the fact that 0;
(k)) 2 > 0.
Therefore, the diagonal matrix & (Es.) + AI 0, substituting which in to (48) yields
II VRx (A(h)) 0.
(49)
Since A is the estimator over F, we have the optimality condition that for any A' E Re x it holds that
max (1 VC (29)
(m)) 11I (120- Al) (1)) s 0.
(50)
Substitute (49) and (50) into P1, for all Go E0 I auc
we have
VC ( (29] (1))
+ +vex (A(i)) + G2GO II, II (12-- )
max A'
[vC ( [2)]as (K)
II, (120-4) A'] (8)
(1m)))
II (VRA (A(a)) II ((20-1) A1] (1))))
max A'
+max
max A'
<0.
(51)
For Term Po.
By definition of V (-), and the regularity condition (v) in Assumption 1, we do the decomposition that U* E V*T, where E0 is diagonal matrix. Projecting V (A(a) into F+ yields that
VQA
(A())
MyI (ve, (Am)) = -Iss. U'UTT) U' (E)) ViT )'T [s. - -V'VIT) -U** --'- (2%.) (v - -V
0.
33

---

High- Siahh-Dimnniinal Tensor Regression with Oracle Properties
Therefore,
+xē") 1151 (12--1)
- AA] (8)) ))
Ps 1x=
Moreover, with the triangle inequality, we have
[ ( (29)
VC (2 ) [[V
1 VL (A4) ((0)
+
VL (1" (A4)
sp
up
s VL
+
[VC
VV (T2(2") (29]
(52)
(16)
where the second inequality comes from the fact that
VC((A)(()) - [vc (
1 [V((A))()) -
vV (2
(R)
sp
From Restricted Smoothness in Assumption 3 where and we have
(AA VC (A + B')IIE - B
E. we can substitute it into (52),
V (A (2))
+L A A 0
VL (AA)))
(53)
(4)
sp
sp
Since 1 (B')= = 0, it is evident that 13 E C. Substitute (42) from Lemma 16 into (53), from the choice of X, we have
!د (IvC (7 (2)1 (m)
VC (29)
sp
(86)
ap
2V WS
V (A* )1())
x (E)
(8))())
sp
n/
ap
S A
By setting = (IVz (A (29)
such that & -UT +2 |
since 20 satisfies the
(k))
nuc
condition 2 E 20 Ftt |20
< 1. we have
1 ([V<(29) (K)
+ +100) =0,
which implies that
P2 = (0 (10 A' Al((x)))= -0.
(54)
Substituting (51) and (54) into (47), we obtain (46) that
([
(k) + 2x (A(a) + +20 [2 A²) (8)
max A'
S0.
Now we are going to prove that 20 A and the error bound between A and A"
Similar to the proof of Lemma 15, since |l-Ilmue is convex, and applying Lemma 13, we have
-A) + (vei (A(s) AG [20- -2 (8) 0 0 VC (A) 22 -
+ (VVC (29) A -=2") + (ve, (Ath) + +220 [a 1A-2") (k)
- -
(p (p-
(55)
F
34

---

High- Dinnnsionnn Tensor Regression with Oracle Properties
From (45), we have
VC (A) A .A-22)) 2²) + (VD, (A()) +2G, (A-2) (13) {(TE(2) A- 1-44)= + vei (2(0))
(21))+AG; +GG, [ A-)) (10) }
S max A'
(56)
Therefore, in (55),
- + (ves (A(x) + +GG, [2 - Al (K) 20.
From (46), we have
(VC (7((2)).2 ..--0 -2) + -2)-) (TO,, (Ais) + +AG. - -A (A)
-
5 max {(Vc(A)) 77 - -A') + (Ve, (At)) + +1G, [2 - -A'] (4))
) 5O.
(57)
A'
Therefore, in (55),
VVC (vc(2)) A 20) +
+ vei (Ai)) + +GG² [ 20) (8)
2 0.
Substituting (55) and (56) into (57) such that
(n- (n- 1|a - 2 2 0.
Since > C, the inequation holds only if
A.
And by Lemma 16, we obtain the statistical oracle bound for the penalty
2//1a
|1> (IX* (8)1
| A-11 1
|
[
A+)
s
(8)
nu
Furthermore, we can have
2v SA
A'
F
np
|1> (Ix* (8))(s))
|gg spp
where Tk-
which completes the proof.
B.5. Proof of Theorem 9
Recall that the proposed alice-wise low- oow-rankness penalty can be reformulated as the sum of the & penalty and concave part. Specifically, we have:
Ummig dm aill
Immagg dmn
R₂ (A) =
Σ
Px (os (IA))))...)))
Σ 1-1
+ ex (1A4)s))...)
1
35

---

High- Diahh-Diinnsinnal Tensor Regression with Oracle Properties
where الهl = min ddydr. Il de). [Aga)].., denotes the 1-th slice of the mode-(j, k) unfolding A(.a) and oo (IAG))))..)) denotes the s-th singular value of the slice. For the estimation problem, we define
Amaaa dm
=(A)=
E
2x (1A(ON))).)))
where Qx (1A6)))),)) = Ea aa (o+ (IA(GA))...))
Based on Lemma 13, for slice-wise lowrankness regularizer, we can similarly prove the following lemmas
||
Lemma 17. Under Assumption 2. p S. and the regularization parameter À
we have
2n
| (12cA0)) ,.s (AGAs) 1..
< $ (1AGSN))... - ) mue
....)
nuc
Proof. Similar the proof of Lemma 14, we can prove the Lemma 17
From Lemma and Lemma 17, we can prove the following general deterministic bound.
Lemma 18. For the estimated parameter tensor À and the true parameter tensor A". we have
1 Umaga dm
||-
|
II
(E))
I.. + +3A S
(
Proof. Similar the proof for Lemma 16, we can derive the error bound for the slice-wise low-rankness regularizer.
Lemma 19. Suppose A" Bdxxxxxxe xdx with rank of each slices |Sp]. Then the error bound between the oracle estimator A and the true A satisfies
I Hmzii dm
2
Se
11=
|[ A
A A'
1
Σ
V
(j.k)
Proof. With Lemma 19, we can also obtain that A0 = A. Similarly, we can prove the Theorem 9.
C. Complementary Experimental Results
In this section, we present additional results for the proposed penalties introduced in Section A. In Section C.I, we evaluate the performance of these penalties on synthetic data and provide a detailed analysis of the experimental findings. Furthermore, Section C.2 demonstrates the effectiveness of the penalties on real-world data.
C.1. Synthetic Data
Figure 5 illustrates the impact of the fiber-wise sparsity regularizer on estimation accuracy. these experiments, we consider 3rd-order tensor A daxdxd, with d 16. We display the results of the Mean Squared Frobenius norm Error (MSFE) when varying the tensor dimension d, the fiber-wise sparsity |Sy]. and the sample size n, respectively. In Figures 5a and 5b, we fix the fiber-wise sparsity |Sal = 4. The three lines in different colors represent varying sample sizes n == (1000, 2000, 3000). From Figure 5a, we observe that MSFE increases as the tensor dimension increases. From Figure 5b, we find that increasing the sample size n decreases the MSFE. This demonstrates that larger sample sizes improve the accuracy of the tensor estimation, as expected. Figure 5c, we see that increasing the fiber-wise-spane wisessparsity |Sg] leads
36

---

High--iimnniionl Tensor Regression with Oracle Properties
0.9
0.45
n1= 1110000
10:000
n-2000
0.8
n--000
0.4
n=0000
0.7
300000
0.35
0.5
0.6
0.3
0.5
0.25
0.4
0.2
0.3
0.15
0.2
0.1
0.1
0.05
10
2000
3000
3500
[Ss]
(a)
(b)
(c)
Figure 5: Fiber-wiss sparsity regularizer with the error bars of MSFE + standard deviation.
0.9
n-1000
n=1000
C4B
n-2000
n~2000
0.4
0.4
n=3000
n-3000
0.6
318PP
MSFE
MSPE
0.2
0.4
0.2
0.3
0.2
0.1
2000
(a)
(c)
Figure 6: Slice-wise sparsity regularizer with the error bars of MSFE + standard deviation
n=-0000
n-1000
d-8
n0:2000
1.5
n=-200
1.2
n=3000
3-0000
d-16
d-20
MSFE
0.8
0.9
10
12
14
15
18
20
2000
2500
5000
3500
4000
|8
(a)
(b)
(c)
Figure 7: Slice-wise low-rankness penalty with the error bars MSFE + standard deviation.
37

---

High-Diimeniooall Tensor Regression with Oracle Properties
Table 2: Comparisons between proposed nonconvex penalties and convex penalties.
Structures
Methods
Synthetic Data
Real-world Data
size
[S]
MSFE
RMSE
MSFE
MPRE
Eitry-wise
Nonconvex
16 x 16 16
2048 0.1
0.1012 + 0.0201
0.0992 + 0.0021
131.5861 # 11.2950
7.6072 + 0.0301
Convex
0.6938 0.0297
0.1004 = 0.0023
141.7100 = 11.9447
7.7198 :: 0.0157
Sparsity
Fiber-wise
Nonconvex
16 16 16
8
0.1
0.4406 +: 0.0157
0.0993 + 0.0012
90.3068 # 7.3006
4.7161 4: 0.0103
Convex
0,7512 + 0.0439
0.0996 + 0.0019
102.7019 + 9:2188
5.0330 + 0.0118
Slico-wisc
Nonconvex
16
8
0.1
0.5761 + 0.0289
0.0997 + 0.0027
43.8706 + 3.0257
1.8909 + 0.0043
Convex
0,7201 + 0.0314
0,1905 + 0.0039
18. 1585 3. 3.8831
1.9250 =: 0.0015
Nanconvex
0.5482 0. 0.0395
0.1002 1 0.0012
35.5536 1: 1.4889
1.0330 1: 0.0022
Mode-wise
16 1616 16
5
Low-rankness
Convex
1.7411
0.0953
0.1096 = 0.0020
41.2719 = 3.5079
1.1027 : 0.0021
Nonconvex
0.9214 + 0.0736
0.1001 + 0.0010
8,9318 + 0.7493
0.0136 + 0.0002
Slice-wise
16 16
x
5 1
Convex
1.8261 0.1056
0.1113 = 0.0031
10.1655 = 0.9050
0.6318 :: 0.0000
to an increase in the estimate error. Furthermore, the standard deviation of the estimation error follows the same trend, increasing with fiber-wise sparsity.
Figure 6 presents the results of the slice-wiss sparsity regularizer. In Figures 6a and 6b, the number of slices is uniformly set to 20. And we set the slice-wise sparsity |Sa] = 4 in Figures 6a and 6c. We select three sample sizes while varying the dimension d or the number of non-zero slices |Sa]] The results indicate that the estimation error increases with increments in d or |Sa]. The standard deviation of the MSFE also rises as the MSFE increases. Furthermore, as observed in Figure 6c, increasing the sample size reduces estimation errors when the dimension is fixed.
Figure demonstrates the results of the Alice-wise low-rankness penalty. In Figure 7b, the r-axis SS|] represents the rank of each slice of the tensor A E pdxdxdd Figures 7a and 7c fix the rank of each slice 5. The three distinct lines correspond to the estimation errors for sample size {1000, 2000, 3000). From Figure 7a, we observe that with a fixed rank and sample size, the estimation error increases as the dimension d enlarges. Furthermore, Figure 7b shows that the estimation errors increase with the rank. Figure 7c demonstrates that with more samples, the estimation errors decrease.
In Table 2, we compare the performance of our proposed nonconvex penalties against traditional convex penalties. For sparsity penalties, we set the r = 0.1, and for low-rankness penalties, we set r) = 1. We configure the tensor dimension such that tensors with wise structures A E Rdxdxx and the others A E dxxxxd where d = 16, 20. Depending on the tensor structure, the sparsity or the rank of the tensors varies accordingly. The results in Table 2 demonstrate that nonconvex penalties achieve lower MSFE for parameter estimation and lower RMSE for predictions compared to their convex counterparts. These empirical findings are in strong agreement with our theoretical analysis.
Table 3: The Frobenius norm
A-
A'
with standard variance changing the noise parameter
Structures
Methods
d d= 10 x 10 x 10
d == 20 x 20 x 20
1 0.1
17 1
17
7) m 0.1
17 1
Entry-wise
1.8909 = 0.2001
1.9021 :: 0.2211
1.9015 + 0.2081
19:1215:2.227
18.6899 + 7556
19.2901 +: 2.3523
Sparsity
Fiber-wise
1.8622 + 0.2813
1.8461 + 0.2783
1.8500 + 0.2431
19.8020 2.6721
19.3542 +: 2.8909
19.7628 + 2.4745
Slice-wise
2,0509 0.3510
2.0121 0.3241
1.9927 0.3666
20.0062 2.7113
20.4267 + 7248
19.4231 + 2.6420
Mode wise
2.6311 0.3008
2.6947
0.3254
2.6265 0.3410
24.6129 2.7010
23.8932
2.7108
24.6571
+: 2,5114
Low-rankness
Slice-wise
3.0150 0.3227
3.0045 4: 0.3754
3.0184 +: 0.3365
25.8502 3.7601
25.6691 3.5732
25.9131 4: 3.8282
Additional experiments. In this paper, we derive five corollaries that establish error bounds involving the noise parameteryy. The analysis of these bounds is nontrivial due to the interplay among conjugate operators, projection operators, and nuclear norm regularization. From equation (20), we observe that increasing 7) enlarges the associated error term, thereby worsening the overall error bound. To illustrate this effect, Table presents results from synthetic data experiments conducted under varying noise levels, which confirm the anticipated impact of 1) on the error magnitude.
As outlined in the five corollaries of our paper, our theoretical framework is inherently generalizable to tensors of any order. Although the scope of this paper did not include experimental results for higher-order tensors, in Table 4, we have conducted supplementary experiments that demonstrate promising outcomes for these cases.
Also, to explore whether the proposed methods perform robustness under an unknown structure, in table 5, we implement
38

---

High- Siahh-Dimnnsinall Tensor Regression with Oracle Properties
Table The Frobenius norm
A
A"
with standard variance for higher dimension F
Struetures
Methods
3-order
4-order
5-order
d
d 16
d 16
d
d 16
Entry-iss
1,0509 1: 0.1004
1.9221 = 0.2214
7.9015 1: 1.2084
19.4215 1: 2:2710
58.6899 1: 7.7566
192.2904 1 17.3523
Spursity
Fiber-wise
1.0622 1. 0.0813
1.8461 0.2783
7.8500 1.1.1.1..4 1.24311
19.8920 1 2.6721
50.3542
1.8. 8.2909
190.7626 1 16.4745
1.0909
1.3066
20,0062
60.4267 8.1248
Mode-wise
1,6311 0.3008
2.6947 0.3254
8.6265
1.3410
34.6129 + 3.4010
63.8932 9.71108
224.6571 + 21.5114
Low-rankness
Slice-wise
1.8150 = 0.3227
2.7045 0.3754
9.0181 + 1.3365
35,8002 = 3.7001
65 6691
+ 9.5732
245.9131 # 22.8282
Table 5: The Frobenius norm |
A+1 with standard variance changing ground data structure of our proposed methods
Structures
Methods
Tensor Data Structures
entry-sp
fiher-sp
slice-sp
Ir-mode
Ir-slice
Entry-wise
1.0509 1: 0.1004
1.0680 1: 0.1027
1.1263 + 0.2046
1.8991 +: 0.4002
1.9367 :: 0.3979
Sparsity
Fiber-wise
1.0931 1 0.1421
1.0622 1 0.0813
1.1305 1 0.2488
1.9054 0.3865
1.9274 1: 0.4410
Slice-wise
1.1014 0.1852
1.1226 = 0.2200
1.0909 0.1510
2.0221
0.4518
2.3185 = 0.4477
Mode-wise
6.8502 1.2101
6.9333 + 1.2565
6.9068 1.1987
1.6311
0.3008
14.9490 = 1.4555
Low-rankness
Slice-wise
7.1481 1 1.3061
7.1636 1 1.2989
7,0701 1 1.3004
15.3770 + 1.3435
1.8150 1 0.3227
Table 6: The average computational time with standard variance comparing nonconvex algorithm and MATLAB CVX solver.
Structures
Methods
8
x 16
Iterations
Total time
Iterations
Total time
Nonconvex (APG Algorithm)
17.5000 0.5415 9.5374
0.2927
27.8345 1 0.3892
16 5965 0.6731 254 251. 1 2.5022
Entry-wise Sparse
Convex (CVX solver)
76.0145 1.2664
Nonconvex (APG Algorithm)
22.3333 + 0.3258
12.5169 0.4709
27.9677 + 0.1196
17.3974 0.5930
Mode-wise Lowrank
Convex (CVX solver)
77.5089 1.0035
229.7236 +: 2.7010
39

---

Biggh-Dimensionll Tensor Regression with Oracle Properties
experiments the proposed methods for each tensor structure, and the results are shown in the table.
C.2. Real-world Data
We have chosen several real-world images from the ImageNet 2012 dataset (Russakovsky et al., 2015) besides the image used in Section 6.2. We implement experiments based on different penalties, revealing the following performance. In Figure& 8, we pick one image "rabbit" from the dataset, and the denoised results are shown in the figure.
(a) "rabbit"
(b) noisy "rabbit"
(c) denoised "rabbit" (cvx. entry-sp)
(d) denoised "rabbit" (nevx. entry-sp)
(e) denoised "rabbit" (cvx. fiber-sp)
(f) denoised "rabbit" (nevx. fiber-sp)
(g) denoised "rabbit" (cvx. slice-sp)
(h) denoised "rabbit"
(nevx. slice-sp)
(i) denoised "rabbit" (cvx. mode-Ir)
(i) denoised "rabbit" (nevx. mode-Ir)
(k) denoised "rabbit" (cvx. slice-Ir)
denoised "rabbit" (nevx. slice-Ir)
Figure 8: The denoising results with the fiber-wise sparsity regularizer.
We have also implemented additional real-world data experiments with the proposed methods. In Table 7, the real data is considered the tensor to be estimated. Regarding the initialization of the covariate tensors A in the real-data experiments, the number of covariate tensors A corresponds to the sample size n = 5000, and the noise term & are drawn independently from a Gaussian distribution with mean O and variance equal to r) = 0.01.
40

---

High- Dithhninnnnsi Tensor Regression with Oracle Properties
Table 7: The MSFE of the climate data (10 years) observation and alcoholic genetic predisposition data with our proposed methods
Sparsity
Low-rankness
Dataset
Penalties
Entry
Fiber
Slice
Mode
Slice
SCAD
7.4556 + 0.8235
8.0716 0.9456
8.4187 + 0.9491
11.6809
2.0437
9.9750 + 1.2203
N4-1990-202-Mommnn Montllllly
MCP
7.6087 + 1.0003
7.9554 + 0.9884
8,1305 + 0.9050
12.6281
2.1882
10.2314 2.0015
Convex
8,2502 + 1.4887
8.7425 # 1.7264
9.4577 + 1.3001
14,5808
3.3435
13.4508 2.6688
SCAD
12.6865 + 2.3544
13.8878 2.2412
14.5610 2.4898
18.7983
5.0917
16.3202 + 4.8331
EEG Dutabuse
MCP
13.0024 + 1.9973
14.0368 + 2.0241
16.0431 + 3. 3.9314
18.4546 + 4.8020
16,4002 4.7771
Convex
13.8001 + 2.6764
14.5716 + 2.1379
16.8890 +: 4.3051
19.6404 5.3317
18.5783 5.3854
The experimental data employed in this study were sourced from the University of Southern California's Viterbi School of Engineering repository and the UCI Machine Learning Repository's 's EEG Database. Specifically, the datasets can be https:
accessed via the following links: https: /ttp://aa..v..aeeeeeeeeeeaaaaeee ics uci //viterbi-web use. www.eeiu3//eeeeeeeeeeeeee data html.
11-1990-2002-Monhhe is monthly climatological dataset (size of 22 x 19500) that includes monthly observations of time series data of 18 climate agents in 125 locations in North America. The original data size for one location in 12 years is a 22 x 156 matrix. To fit our model, we segment it into twelve 22 13 matrices and merge them into 22 x 13 x 12 tensor to predict. The estimation results are shown in Table 7.
EEG Database is an alcoholic genetic predisposition dataset that contains EEG images of 64 channels sampled at 256 Hz for subjects suffering from alcoholism and 44 normal controls. In the dataset, there are 10 alcoholic subjects, and each sample is a third-order tensor (Channelsx Time Voltage). We take each sample the tensor estimate, and the result revealed in Table 7.


Reviews:
[{"review_id": "3Rvkl1UU3J", "reviewer": "Reviewer_xpr9", "summary": "This paper proposes a framework for tensor on tensor regression, introducing novel nonconvex regularizers and an accelerated proximal gradient algorithm for estimation. The authors propose a class of penalties that depend on the singular values of each tensor dimension and give Frobenius norm rates of convergence guarantees under oracle optimal hyper parameter tuning. A proximal gradient algorithm is also provided as a feasible estimation procedure. Numerical experiments and an empirical application show the advantage on the proposed method.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely appreciate your positive comments and insightful feedback. We have carefully addressed your concerns as follows:\n\n**1. Justification for Claim 2**\n\nIn this paper, we focus on nonconvex sparse learning in tensor regression. Nonconvex regularizers such as SCAD and MCP can achieve the oracle estimation rate because they induce less bias compared to convex penalties like the $\\ell_1$ norm. Unlike Lasso, which uniformly shrinks all coefficients and can introduce significant bias for large signals, nonconvex penalties apply little to no shrinkage to large coefficients while still promoting sparsity among small ones. This selective regularization enables accurate support recovery and nearly unbiased estimation on the true support, leading to improved statistical efficiency under appropriate conditions.\n\n\nIn Corollary 1, we present the estimation performance for the sparse parameter under nonconvex regularization. We also compare this result with the estimation rate under convex regularization. The corollary explicitly shows that nonconvex regularization enables improved estimation accuracy by eliminating dependence on the ambient dimension $d$. For the low-rank parameter estimation setting, the corresponding result is provided in Corollary 2, which can be directly compared with prior work such as [1,2]. For instance, Lemma 10 of [1], which uses convex nuclear norm regularization, yields an estimation bound that scales with the tensor dimension. In contrast, our bound remains dimension-independent.\n \n\n[1] G. Raskutti, M. Yuan, and H. Chen. Convex regularization for high-dimensional multiresponse tensor regression. In Proceedings of the 36th International Conference on Machine Learning (ICML), 2019.\n\n[2] S. Negahban and M. J. Wainwright. Estimation of (near) low-rank matrices with noise and high-dimensional scaling. Annals of Statistics, 39(2):1069–1097, 2011.\n\n\n**2. Comparison with Rate from LASSO**\n\nThe estimation error under nonconvex regularization differs from—and indeed improves upon—that of Lasso, as it achieves a faster convergence rate equivalent to the rate attainable when the true support is known in advance. Specifically, by leveraging nonconvex penalties, the estimation bound becomes independent of the ambient dimension $d$, thereby achieving the oracle rate. This dimension-free behavior highlights a key advantage of nonconvex regularization.\n\n\n\n**3. When Might Convex Methods Outperform Nonconvex Methods?**\n\nThank you for this question. While nonconvex regularization can yield improved estimation rates under ideal conditions, convex methods such as Lasso may outperform them when these assumptions are not satisfied. For instance, nonconvex methods often rely on strong signal conditions to ensure accurate support recovery. In low signal-to-noise ratio settings, where such conditions may fail, convex methods can produce more stable and reliable estimates due to their uniform shrinkage behavior and greater algorithmic robustness.\n\n\n**4. Clarification on Table 1** \n\nDue to space constraints, we report results for two representative nonconvex regularizers in Table 1. Results for additional regularizers are provided in the Appendix. We will include more detailed discussions of both Table 1 and Table 2 in the camera-ready version.\n\n\n**5. On Theorem 5** \n\nTheorem 5 provides a general characterization of the statistical estimation performance under nonconvex regularization in tensor regression. Notably, the choice of regularizer—convex or nonconvex—affects both the Gaussian width and the subspace compatibility constant, thereby highlighting fundamental differences between the two approaches. For specific penalty choices, the subsequent corollaries clearly demonstrate that nonconvex regularization can achieve faster convergence rates than convex counterparts.\n\nBesides, we will consider the suggestion to rename Theorem 5. \n\nWe hope the above responses address your comments clearly, and we sincerely thank you again for your valuable feedback and thoughtful review of our paper.", "claims_and_evidence": "The paper makes the following claims:\n\n1. The proposed nonconvex penalty estimators for tensor on tensor regression achieve oracle optimal convergence rates under different sparsity assumptions. The estimator can be computed through a proximal gradient algorithm.\n\n2. The nonconvex penalty estimators exhibit faster convergence rates compared to those with convex penalties.\n\nOverall the paper makes a very strong case for claim 1. Claim 2, to my reading of the paper, is mostly supported by the empirical exercises as the theoretical results are not directly compared to the convex penalty cases. With this mind, expanding on the justification for claim 2 would improve the value added of the paper.", "methods_and_evaluation_criteria": "The methods and evaluation criteria are adequate for the problem at hand.", "theoretical_claims": "The theoretical claims are well stated and appear novel and correct.\n\n* Given the claims of the paper it would be nice if more direct theoretical comparisons would be made with convex penalty estimators. Is there a theoretical result that ensures that the rates of convergence are faster for the nonconvex penalty? It seems from corollary 1 and 2 that the rates are very similar to the cases for LASSO for example.", "experimental_designs_or_analyses": "The analysis are careful and convincing. However, more detail on the comparison with convex methods in table 1 and the simulation exercise would be nice to reinforce the claims of the paper. In which cases do we expect the convex methods to perform better? Expanding the simulation exercises to study this would improve the relative contribution of the paper.", "relation_to_broader_scientific_literature": "I am not very familiar with the tensor regression literature but it would be good to clarify the novelty of the results and Theorem 5 vis a vis the literature, in particular to clarify why the rates are faster than for the convex penalty estimators.", "essential_references_not_discussed": "-", "other_strengths_and_weaknesses": "Overall the paper is very well written and makes a very compelling case for their proposed method! My only concern is whether the authors actually provide a justification for the nonconvex penalties having a faster convergence rate, beyond the empirical simulations.", "other_comments_or_suggestions": "* Theorem 5 should be renamed Theorem 1.", "questions_for_authors": "See above.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "V3vQ1xBLgG", "reviewer": "Reviewer_8rnB", "summary": "This paper studies tensor regression models with non-convex penalties and provides general risk bounds for the resulting estimators, demonstrating that they achieve oracle statistical rates under mild technical conditions. The authors also introduce an accelerated proximal gradient algorithm to estimate the proposed estimators.  \n\nA comprehensive set of experiments is conducted to illustrate the advantages of non-convex penalties over convex penalties. While the theoretical results are solid, the paper primarily extends [1] to the non-convex penalties. As a result, the novelty is somewhat limited, and certain aspects of the presentation are unclear.  \n\n[1] Raskutti, G., Yuan, M., & Chen, H. (2019). Convex regularization for high-dimensional multiresponse tensor regression.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**More discussions on decomposition-based tensor regression methods**\n\nWe appreciate the comment regarding tensor regression methods based on tensor decomposition techniques such as CP and Tucker. Both decomposition-based and regularization-based approaches have been actively explored in the tensor regression literature, and we have reviewed these works in the Introduction. However, we would like to emphasize that the primary focus of this paper is theoretical—namely, to establish oracle statistical guarantees for tensor regression under nonconvex regularization. To the best of our knowledge, existing decomposition-based methods have not addressed this theoretical direction. Therefore, decomposition-based approaches represent a distinct line of research, and direct comparisons fall outside the intended scope of this work.\n \n**Experiment on Real-world Data**\n\nIn our experiment, we did not have access to a dataset with $\\mathcal{X}$ and $\\mathcal{Y}$ as direct measurements. Therefore, we generated such a dataset based on a ground truth tensor image $\\mathcal{A}^\\ast$, which is a common practice in the literature (see [1,2]). This setting allows us to directly evaluate the estimation performance by comparing the estimated tensor $\\hat{\\mathcal{A}}$ with the ground truth $\\mathcal{A}^\\ast$. Additionally, we have also reported the error between $\\mathcal{Y}$ and $\\widehat{\\mathcal{Y}}$, which assesses the model's predictive performance when $\\mathcal{A}^\\ast$ is unknown—precisely the scenario raised in your comment.\n\n[1] Romera-Paredes, B., H. Aung. “Multilinear multitask learning”. In International Conference on Machine Learning.\n\n[2] Liu, Y., Liu, J., Long, Z. Tensor regression. Springer.\n\n**Response to weaknesses**\n\n1. In this paper, the term *oracle* refers to the performance of an estimator that assumes knowledge of the true support, as in [3,4]. For example, in the case of element-wise sparsity, the oracle estimator $\\widehat{\\mathcal{A}}^{O}$ satisfies $\\Vert\\widehat{\\mathcal{A}}^{O}-\\mathcal{A}^*\\Vert_{\\mathrm{F}}\\lesssim\\Vert(\\nabla L(\\mathcal{A}^*))_ {\\mathcal{S}^*}\\Vert_{\\mathrm{F}} \\asymp \\sqrt{\\frac{s}{n}}$, which follows directly from a first-order Taylor expansion via the mean value theorem.\n\n[3] Fan,J.,Liu,H.,Sun,Q. I-LAMM for sparse learning: Simultaneous control of algorithmic complexity and statistical error. The Annals of Statistics.\n\n[4] Gui,H.,Han,J. Towards faster rates and oracle property for low-rank matrix estimation. In International Conference on Machine Learning.\n\n2. We will correct the reference issue in the camera-ready version.\n\n3. [Explanation on Figure 3] Figure 3 visualizes the element-wise estimation error between a randomly generated third-order tensor and its estimate obtained using our nonconvex low-rank regularizer. Each element is evaluated individually: if the absolute error exceeds a fixed threshold, the point is marked in red; otherwise, it is shown in blue. Figures 3(a) and 3(d) present representative outcomes using the nonconvex method, while Figures 3(b) and 3(e) illustrate results from a convex regularization approach. We observe that the nonconvex method produces significantly fewer red points, indicating smaller element-wise errors and demonstrating the improved accuracy achieved by our proposed approach.\n\n**Response to Questions**\n\n1. All reported experimental results are based on 100 Monte Carlo replications, as stated in the first paragraph of Section 5. We will update this to 1,000 replications in the final version to ensure greater robustness and stability of the results.\n\n2 & 3. We acknowledge the reviewer’s concern regarding the anomalous behavior where the standard deviation appears to increase with the sample size. We sincerely appreciate this observation. After a careful review, we found that this issue was due to an insufficient number of Monte Carlo replications. To address this, we conducted additional experiments with 1,000 Monte Carlo trials, which yielded more stable and robust results. The updated figures, including standard deviations and computational time measurements, are available at the following link: [https://anonymous.4open.science/r/ICML_Rebuttal-6182/](https://anonymous.4open.science/r/ICML_Rebuttal-6182/). We thank the reviewer again for this helpful comment and will correct the issue in the final version.\n\n**Essential References Not Discussed**\n\nWe are grateful to the reviewer for bringing this valuable reference to our attention, and we will cite it appropriately in the final version of the paper.\n\nThank you for your thoughtful and constructive feedback. We sincerely appreciate your advice and the time you took to review our work. However, we are somewhat puzzled by the decision to reject the paper. We hope that the above clarifications and updated results have addressed your concerns. If our responses have sufficiently resolved the issues raised, we would be grateful if you could consider reflecting this in your final evaluation.", "claims_and_evidence": "The authors present the primary theoretical results for the general estimator in Equation (2) and derive upper bounds under different scenarios. However, I find the explanation of the oracle property unclear. In Section 3.2.1 SPARSITY REGULARIZATION, the authors first define the oracle rate as the statistical convergence rate of the oracle estimator and provide a detailed element-wise sparse oracle estimator, which is well-defined. However, in Corollary 1, the authors present the rate of $\\hat{\\mathcal A}$ without establishing any explicit connection to the defined oracle estimator $\\widehat {\\mathcal A}^O$. It is unclear how the authors conclude that $\\hat{\\mathcal A}$ achieves the oracle rate under some weak assumptions. A similar issue arises in Section 3.2.2 **LOW-RANK REGULARIZATION**.  \n\nTo rigorously establish the oracle property, the authors should provide theoretical results like that:  \n1. $\\hat{\\mathcal A}_{\\bar S_1}=0$\n2. the asymptotic normality of ${\\hat{\\mathcal A}_{S_1}}$ \n\nAdditionally, the paper includes extensive experiments to demonstrate the advantages of non-convex penalties over convex penalties. However, in the real data experiment, the authors assume that the ground-truth $\\mathcal A$ is known and simulate $\\mathcal X$ and $\\mathcal Y$ accordingly. A more realistic evaluation would involve an experiment where only the observed $\\mathcal X$ and $\\mathcal Y$ are available, allowing an assessment of the model's performance when $\\mathcal A$ is unknown.", "methods_and_evaluation_criteria": "In the experimental section, the authors primarily compare convex penalties. However, there exist many other comparative methods for tensor regression on a single dataset, such as tensor regression based on CP decomposition and Tucker decomposition.", "theoretical_claims": "I have read the proof in the author's appendix, which is very detailed. I believe it is feasible, but the specific details require further reading.", "experimental_designs_or_analyses": "I believe the soundness/validity of any experimental designs or analyses is feasible. The authors have conducted experiments on both simulated and real data. However, the experiments provided so far are all based on simulations, assuming the true tensor coefficients are known. The authors should include real data analysis using actual tensor covariates and response \n for analysis (without knowing the true ). This would make the findings more convincing.", "relation_to_broader_scientific_literature": "The theoretical analysis primarily compares the proposed method with convex penalties, particularly emphasizing that the use of non-convex penalties achieves the oracle rate. However, the current presentation lacks clarity in this aspect.", "essential_references_not_discussed": "The authors could further discuss more papers about regularized tensor regression based on tensor decomposition methods, for example, [2].\n\n[2] Lu, W., Zhu, Z., & Lian, H. (2020). High-dimensional quantile tensor regression. Journal of Machine Learning Research, 21(250), 1-31.", "other_strengths_and_weaknesses": "### Strengths  \nAs discussed above.  \n\n### Weaknesses  \n1. The discussion on the oracle property is unclear. It is difficult to directly see how $\\hat{\\mathcal A}$ and the oracle estimator exhibit similar performance.  \n2. There are issues with the citations. For example, *(Hua Zhou & Zhu, 2013)* and *(Zhou et al., 2013)* refer to the same paper. The authors should carefully check and correct the references.  \n3. The authors should report the computational time in the experiments. Additionally, more detailed explanations should be provided for the tables and figures. For instance, **Figure 3** is difficult to interpret intuitively.", "other_comments_or_suggestions": "The authors could consider reorganizing the structure of the paper to better highlight its key contributions. The most significant contribution of this work lies in its theoretical advancements. Therefore, the authors should emphasize this aspect more prominently and provide a clearer discussion of the challenges in the theoretical analysis and how they are addressed.", "questions_for_authors": "1. In **NUMERICAL EXPERIMENTS**, how many runs were performed for each setting?  \n2. In the numerical results, the standard deviation often increases as the sample size grows. For instance, in **Figure 1a**, when $d=20$, the standard deviation is largest at $n=3000$. A similar trend is observed in **Figure 2b**. Is this due to an insufficient number of runs?  \n3. Can the authors provide a comparison of the computational time across different methods?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "rgOgukWAGA", "reviewer": "Reviewer_1JFB", "summary": "This paper addresses tensor regression models for high-dimensional tensor data. Specifically, it proposes a tensor-response tensor regression model, assuming low-dimensional structures such as sparsity and low rankness. While conventional convex penalties are easy to optimize, they often fail to model the data accurately. On the other hand, non-convex penalties have a higher potential to model the data correctly, but there's no guarantee of reaching the optimal solution. To address this issue, the authors theoretically prove that non-convex regularization terms, such as SCAD (Smoothly Clipped Absolute Deviation) and MCP (Minimax Concave Penalty), exhibit oracle performance under certain assumptions. They then propose an algorithm based on the accelerated proximal gradient algorithm to efficiently compute estimators based on these non-convex penalties. Evaluation experiments using synthetic and real-world datasets demonstrate the effectiveness of the proposed regression model and the practicality of the theoretical results.\n\n## update after rebuttal\n\nThank you for your response. I understand now. I believe the contribution of this paper is significant, so I will keep my evaluation as it is. I mistakenly posted my comment in the Official Comment section. My apologies.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your positive feedback. \n\n\n**Hyperparameters for the Convex Penalties**\n\nAll hyperparameters are selected via ten-fold cross-validation, as noted in the first paragraph of Section 5. This includes the regularization parameter $\\lambda$ used in convex penalties. Specifically, $\\lambda$ is chosen from a uniform grid of 21 values in the range $[0, 1]$, and the optimal value is selected through ten-fold cross-validation.\n\n\n**Other Comments or Suggestion 1 and 2**\n\nWe will correct typos in the camera-ready version.\n\n\nWe greatly appreciate your recognition of our work. If our responses have sufficiently addressed your concerns, we would be grateful if you could consider reflecting this in your evaluation. Thank you once again for your time and valuable feedback.", "claims_and_evidence": "It is intuitively understandable that convex penalties cannot effectively capture data characteristics, especially in high-dimensional tensor data. The authors rigorously prove, through sound theoretical development, that estimators based on certain non-convex penalties exhibit oracle performance under specific conditions. Furthermore, they demonstrate, using both synthetic and real-world data, that non-convex penalties can indeed effectively capture the underlying features of the data.", "methods_and_evaluation_criteria": "Regarding this paper's key claim, Oracle performance, the evaluation appropriately assesses performance variations under different knowledge conditions. While a comparative evaluation against convex penalties is presented, the method for determining hyperparameters when using convex penalties is not explicitly stated. It is presumed that a 10-fold cross-validation approach was employed, similar to the proposed method. Explicitly stating this would enhance the persuasiveness and rigor of the evaluation.", "theoretical_claims": "I have reviewed the proofs, particularly those demonstrating that the non-convex penalty achieves Oracle performance under specific conditions, and found the theoretical development sound and without issue.", "experimental_designs_or_analyses": "As mentioned in the \"Methods and Evaluation Criteria\" section, the evaluation appropriately assesses performance variations under different knowledge conditions, which is directly relevant to the paper's central claim of oracle performance. The evaluation metrics used are also deemed suitable.", "relation_to_broader_scientific_literature": "While tensor completion with low-rank constraints is a widely studied and effective technique, this paper makes a significant contribution to the field by focusing on the potential limitations of convex penalties and providing theoretical guarantees for estimators based on non-convex penalties, which offer the potential for more appropriate modeling.", "essential_references_not_discussed": "The paper appropriately cites relevant prior work on various constraints used in tensor regression. I do not have any essential references to suggest for addition.", "other_strengths_and_weaknesses": "As previously mentioned, this paper makes a significant contribution to the field of tensor learning with various constraints, such as low-rank constraints, by providing theoretical guarantees for estimators based on non-convex penalties. The evaluation experiments are generally comprehensive. However, a concern remains regarding whether the hyperparameters for the convex penalties, used as comparison methods, were optimally determined.", "other_comments_or_suggestions": "There are a few potential errors in the manuscript:\n\n1. In Section 2.1, while $\\mathcal{A}^{\\*}$ is defined as an M-th order tensor and $\\mathcal{X}$ as an N-th order tensor, the expression around line 85 appears to represent the elements of $\\mathcal{A}^{\\*}$ with N indices and the elements of $\\mathcal{X}$ with M indices. It seems that $\\mathcal{A}^{\\*}$ and $\\mathcal{X}$ may be incorrectly assigned in this expression.\n\n2. In the first sentence of Assumption 4, $X^{n)}$ appears. Should this be corrected to $X^{(n)}$, adding the missing left parenthesis?", "questions_for_authors": "How were the hyperparameters for the convex penalties, used as comparison methods in Table 1, determined?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "SOCZ2Zvgoy", "reviewer": "Reviewer_k4E1", "summary": "The paper introduces a high-dimensional tensor-response tensor regression model under low-dimensional structural assumptions, such as sparsity and low-rankness. \n\nThe authors propose a least squares estimation framework with non-convex penalties and derive general risk bounds for the resulting estimators. \n\nThe paper also particularizes the bounds for the case where the support of the solution is known (oracle estimator). \n\nAn accelerated Proximal Gradient Algorithm is proposed and tested on synthetic data and an image denoising problem.\n\n## update after rebuttal\nI thank very much the authors for addressing my questions. However, I keep my original score.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Broader Scientific Literature**\n\nThank you for introducing these relevant references. They are very helpful, and we will cite them appropriately in the final version.\n\n**Response to Claims and Evidence**\n\n1. To the best of our knowledge, the proposed bound is tight in terms of its order, with certain constant factors omitted. Since our analysis focuses on the convergence rate rather than explicit constants, directly comparing the theoretical bounds with observed experimental errors can be challenging.\n\n2. The choice of regularization penalties affects the bound primarily through their impact on the Gaussian width and the subspace compatibility constant, as outlined in Theorem 5. Detailed results for both sparse and low-rank regularization are presented in the subsequent corollaries. For example, in the case of element-wise sparsity, the key difference between convex and nonconvex penalties is that nonconvex regularization leads to improved error rates by reducing the dependence on the ambient dimension $d$, thereby attaining the oracle rate of $\\sqrt{\\frac{s}{n}}$.\n\n3. Assumption 4 on sub-Gaussianity is a standard condition widely used in high-dimensional statistical estimation (see [1]). It is a mild assumption that holds in many practical scenarios and covers a broad class of commonly encountered distributions. We acknowledge, however, that there are settings beyond the sub-Gaussian framework—for example, when the data exhibit heavy-tailed behavior. Extending the theoretical analysis to such cases is an important direction in statistical estimation and can be left for future research.\n\n[1] Raskutti, G., Yuan, M., & Chen, H. (2019). Convex regularization for high-dimensional multiresponse tensor regression.\n\n**Response to Weaknesses**\n\n1. Due to space constraints, only a subset of simulation results is reported in the main text. In fact, many additional experimental results are provided in the appendix, offering a more comprehensive analysis.\n\n2 & 3. Please refer to the reply above.\n\n4. It is well acknowledged in the statistical literature that nonconvex penalties can yield better estimation performance than convex ones by reducing the bias introduced by regularizers such as the $\\ell_1$-norm. In this paper, we investigate how such improvements are reflected in the context of tensor regression. Specifically, nonconvex penalties impose less shrinkage on large coefficients, enabling more accurate recovery of the true signal. This often results in faster convergence rates, particularly under sparsity or low-rank assumptions, and can achieve the oracle rate under suitable conditions.\n\n**Additional Experimental Evaluations.**\n\nThank you for this question. Neuroimaging is indeed an important application area for tensor regression. In fact, we have already conducted relevant experiments. Due to space constraints, the results are provided in the appendix. Our experiments include an electroencephalography (EEG) dataset, which is a form of neuroimaging [2]. We appreciate your suggestion and plan to apply our methods to additional neuroimaging datasets in future work.\n\n[2] Liu, Y., Liu, J., Long, Z., and Zhu, C. Tensor regression. Springer, 2022.\n\nWe sincerely appreciate your recognition of our work. If our responses have adequately addressed your concerns, we would be grateful if you could consider reflecting this in your evaluation. Thank you once again for your time and thoughtful feedback.", "claims_and_evidence": "The paper explicitly or implicitly makes the following key claims:\n\nTheoretical Contributions: The main theoretical result is Theorem 5, which provides an upper bound on the Frobenius norm of the estimated tensor A. While this result is valuable, its practical implications are not thoroughly explored. For instance:\n\n-\tIs the bound tight? In other words, how does it compare with the actual errors observed in the experiments?\n\n-\tHow does the nonconvexity of the penalty affect this bound? Can we compare the bounds for the convex and nonconvex cases?\n\n-\tAssumption 4, which assumes multivariate normality in the data, may not be realistic for real-world applications.\n\nProposed Framework: The paper introduces a general framework for linear regression with tensor-structured input and output data, extending beyond traditional algorithms that primarily handle scalar, vector, or matrix data. This theoretical framework is well-constructed and allows for the incorporation of low-dimensional priors, such as low rank and sparsity. However, the experimental validation is limited to: A simple synthetic dataset with scalar outputs; and a basic denoising problem involving real images.\n\nNonconvex vs. Convex Penalties: The results indicate that nonconvex penalty functions outperform convex ones, as suggested in prior literature. The authors demonstrate that their method, when using a nonconvex penalty, achieves lower error compared to a classical convex penalty based on the ℓ1-norm. However, the paper lacks a theoretical explanation or intuition as to why nonconvex penalties yield better performance.", "methods_and_evaluation_criteria": "The theoretical methodology is sound. The experimental methodology used to evaluate the proposed algorithms is valid but somewhat limited (see detailed comments below).", "theoretical_claims": "The proof of the theoretical claims are provided in the supplementary material, which I did not check carefully. However, the theoretical claims sounds reasonable.", "experimental_designs_or_analyses": "Yes, I check the experimental results and analysis provided within the paper and the supplementary material. The proposed experiments are correct but limited (see detailed comments below).", "relation_to_broader_scientific_literature": "Although the proposed algorithm is evaluated only in limited scenarios—specifically on synthetic datasets and a simple denoising problem—it has the potential to drive significant advancements in various scientific fields where tensor-structured datasets are prevalent. In particular, it could be highly valuable in domains such as neuroimaging, where modeling relationships between input and output tensor data is crucial, as well as in other fields that rely on tensor-based linear regression. Expanding the evaluation to more diverse and realistic datasets would further strengthen the impact and applicability of the proposed approach.", "essential_references_not_discussed": "The paper includes relevant previous references on tensor regression. However, there are important previous works addressing similar problems as it is the case of the following papers: \n\n-\t\"Higher-Order Partial Least Squares (HOPLS): A Generalized Multi-Linear Regression Method\", Q Zhao, CF Caiafa, DP Mandic, ZC Chao, Y Nagasaka, N Fujii, L Zhang, A Cichocki, IEEE Trans. on Pattern Analysis and Machine Intelligence (PAMI), Vol. 35 No. 7, 2013. doi:/10.1109/TPAMI.2012.254.\n\n-\t“A Multilinear Subspace Regression Method Using Orthogonal Tensors Decompositions”, Q Zhao, CF Caiafa, DP Mandic, L. Zhang, T. Ball, A. Schulze-Bonhage, A. Cichocki, Proc. NIPS 2011 (Neural Information Processing Systems), Granada, Spain, 12-17 December 2011.", "other_strengths_and_weaknesses": "Strengths:\n\n-\tThe paper presents a general theoretical framework for linear regression with input and output tensors of arbitrary dimensions, extending beyond traditional approaches that focus on scalar, vector, or matrix data.\n\n-\tA theoretical bound on the estimator error is derived, providing a rigorous foundation for analyzing the performance of the proposed method.\n\n-\tThe study explores the use of nonconvex penalty functions, which appear to play a crucial role in enhancing performance, particularly for high-dimensional datasets.\n\nWeaknesses:\n\n-\tLimited experimental validation: The evaluation is restricted to a simple synthetic dataset with scalar outputs and a basic denoising problem using real images, limiting the demonstration of the framework’s broader applicability.\n\n-\tLack of practical insights from Theorem 5: The paper does not thoroughly explore the practical implications of the theoretical bound. For example, is the bound tight? How does it compare with actual errors observed in experiments?\n\n-\tUnrealistic data assumption: Assumption 4 relies on multivariate normality, which may not hold in real-world applications, potentially affecting the method’s generalizability.\n\n-\tMissing theoretical justification for nonconvex penalties: While experiments show that nonconvex penalties outperform convex ones, the paper lacks a theoretical explanation or intuition behind this behavior.", "other_comments_or_suggestions": "I have not found typos", "questions_for_authors": "Regarding Theorem 5:\n\n- Is this bound tight? In other words, how this bound is compared with the actual errors attained in the experiments?\n\n\n- What is the effect of the non-convexity of the penalty on this bound? Can we compare the bounds for the convex and nonconvex cases?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: The reviewers find the paper well written and the claims well supported by experimental evidence. They further find the methodology sound and the theoretical contribution that certain non-convex regularizations in the context of tensor regression exhibit oracle properties a strongpoint.  The reviewers further find that the contribution can provide significant advancements in fields relying on tensor-structured datasets. 

Whereas the initial material in regards to the experimentation and theoretical benefits of non-convex vs. convex regularizations was found to have room for improvements, the authors well addressed the reviewers concerns. In particular, whereas the oracle property was found unclear the authors very well addressed this concern providing in the rebuttal also a proof sketch of strong oracle properties that they promise to include in a camera ready version of the manuscript. During the discussions, this was highly appreciated and deemed an impressive additional theoretical contribution to the paper strengthening the contribution. 

The reviewers agreed the manuscript warrants publication at ICML and the authors are strongly encouraged to carefully update their camera ready version with the substantial improvements and additional clarifications, experimentations and new proof of a strong oracle property provided in their rebuttals in the final version of the paper.

(Note also when correcting minor typos for the final version the following very minor additional typo: “There exits -> There exists”)
Decision: Accept (poster)