arXiv CSopen access

Guaranteed Low-Rank Tensor Recovery from Modewise Measurements via Normalized Block-Weighted Riemannian Gradient Descent

· arxiv_cs

arXiv CS · Papers · License: Open Access

machine learning, deep learning, neural networks

Guaranteed Low-Multilinear-Rank Tensor Recovery from Modewise Measurements via Normalized Block-Weighted Riemannian Gradient Descent * Yushi Zhou†

Feng Zhang‡

arXiv:2609.24679v1 [cs.LG] 21 Sep 2026

September 22, 2026

Abstract We consider the recovery of low-multilinear-rank tensors from linear measurements and propose an adaptive block-weighted modewise Riemannian gradient descent method. The method combines memory-efficient modewise measurements with a normalized adaptive weighting strategy for the core and factor components of the Riemannian gradient. The weighting improves convergence without increasing the multilinear-rank bound of the search direction or the size of the reduced core used for retraction. Under the tensor restricted isometry property and a suitable initialization, we establish local linear convergence and derive sampling guarantees for sub-Gaussian and subsampled orthogonal with random sign (SORS) measurements. Numerical experiments on synthetic low-Tuckerrank tensors show that the proposed method reduces iteration counts and computational time while maintaining reliable recovery performance, especially near the recovery threshold and for structured SORS measurements.

Introduction

The tensor recovery problem arises in a wide range of applications, including computer vision [1, 2], image inpainting [3, 4], multi-temporal image reconstruction [5], machine learning [6], signal processing [7, 8], and quantum state tomography [9, 10]. In this paper, we consider the recovery of an unknown tensor X ∈ Rn1 ×···×nd that admits a Tucker decomposition [11] of multilinear rank r = (r1 , . . . , rd ), from a limited number of linear measurements y = L (X ),

(1)

Q where L : Rn1 ×···×nd → RM is a linear measurement operator and M ≪ di=1 ni . Since the measurement dimension is considerably smaller than the ambient tensor dimension, the recovery problem is generally underdetermined without additional structural assumptions. By imposing the low-multilinearrank prior, tensor recovery can be formulated as the constrained least-squares problem min

X ∈Rn1 ×···×nd

1 ∥L (X ) − y∥22 2

subject to

mulrank(X ) = r.

(2)

A conventional realization of the measurement operator Qd L in (1) first vectorizes the tensor and then applies a single dense measurement matrix of size M × i=1 ni . This vectorization-based construction * This

work was supported in part by the National Key Research and Development Program of China under Grant 2023YFA1008502; in part by the Fundamental Research Funds for the Central Universities under Grant SWU-KR25013; and in part by the National Natural Science Foundation of China under Grant 12101512. † School of Mathematics and Statistics, Southwest University, Chongqing, China. Email: [email protected]. ‡ Corresponding author. School of Mathematics and Statistics, Southwest University, Chongqing, China. Email: [email protected].

underlies many classical algorithms for solving the recovery problem (2) and its sparse counterpart, such as ℓ1 -minimization [12, 13], CoSaMP [14, 12], and iterative hard thresholding [15]. Its practical drawback, however, lies in the size of this matrix: for high-order or large-scale tensors, merely generating and storing this matrix can require more memory than storing the target tensor itself. Consequently, although dense vectorized measurements provide a convenient theoretical model, their storage cost can render them impractical for large-scale tensor recovery. Modewise measurement operators provide a structured and memory-efficient alternative to conventional vectorization-based measurements. Following the modewise measurement framework of [16, 17], as further developed in [18], a two-stage modewise linear operator L : Rn1 ×···×nd −→ Rm1 ×···×md′ can be expressed as L (X ) := R 2 R 1 (X ) ×1 A1 ×2 · · · ×de Ade ×1 B1 ×2 · · · ×d′ Bd′ , (3) where R 1 : Rn1 ×···×nd −→ Rne1 ×···×ende is a reshaping operator that reorganizes the entries of the original e i ×e ni , i ∈ [d], e e is applied along d-mode tensor into a d-mode tensor. At the first stage, each matrix Ai ∈ Rm m e 1 ×···×m e de the corresponding mode of the reshaped tensor, yielding an intermediate tensor in R . The opm′1 ×···×m′d′ m e 1 ×···×m e de erator R 2 : R −→ R subsequently reorganizes the first-stage output into a d′ -mode ′ tensor. At the second stage, the matrices Bj ∈ Rmj ×mj , j ∈ [d′ ], are applied along the corresponding modes to further compress the tensor, producing the final measurement tensor in Rm1 ×···×md′ . This framework can be extended to construct more general multi-stage operators by repeatedly alternating the reshaping and modewise-compression procedures. Conventional vectorized measurements can be viewed as a particular case of (3), in which R 1 vectorizesQthe entire tensor and the first stage compression is represented by a single dense matrix of size M × di=1 ni . In contrast, when R 1 performs only a moderate reshaping, the modewise construction replaces this massive matrix with a collection of substantially smaller component matrices. It therefore requires fewer random parameters and reduces the storage cost of the measurement operator. The individual mode products are also naturally amenable to parallel implementation and preserve the multimodal organization of the tensor throughout the first stage compression. The storage savings and structural advantages noted above are also characteristic of related approaches, including modewise oblivious subspace embeddings, Kronecker-structured Johnson–Lindenstrauss transforms, and sketching methods for large-scale Tucker approximation [16, 19, 20]. A key theoretical ingredient that connects such structured embeddings with low-rank tensor recovery is the tensor restricted isometry property (TRIP). Introduced in the analysis of tensor iterative hard thresholding, TRIP extends the classical restricted isometry principle to low-rank tensor models by requiring the measurement operator to approximately preserve the Frobenius norm of tensors with prescribed rank structure [21]. It therefore provides a natural condition under which the geometry of the low-rank tensor model is retained after measurement. In the modewise setting, when the component matrices satisfy suitable restricted isometry conditions, the resulting one-stage and two-stage operators likewise satisfy TRIP over tensors with bounded multilinear rank [18]. In particular, these guarantees can be established when the component matrices are chosen from either subgaussian ensembles or Subsampled Orthogonal with Random Sign (SORS) ensembles, the two measurement constructions considered later in this work. Motivated by these results, we specialize the general framework (3) to the one-stage operator L 1 = vec ◦A ◦ R and the two-stage operator L 2 = A2nd ◦ L 1 , and incorporate them into the proposed Riemannian recovery method. The corresponding TRIP-based convergence and sampling guarantees are established in Section 4. Related work A variety of methods have been developed for low-rank tensor recovery. Convex approaches commonly promote low multilinear rank by minimizing the sum of the nuclear norms of tensor matricizations [22, 23, 24], but matricization may fail to fully exploit the intrinsic multilinear structure of the tensor and can lead to suboptimal sampling complexities. Tensor nuclear-norm formulations provide a more direct alternative [25], but evaluating or optimizing the corresponding tensor norms is generally computationally demanding. More recently, nonlocal tensor nuclear norms have been developed to 2

exploit spatial correlations in high-dimensional image recovery [26]. Nonconvex tensor low-rank surrogates have also been employed in imaging inverse problems, such as tensor logarithmic Schatten-p minimization for 3D Poissonian image deblurring [27]. Nonconvex approaches have therefore received considerable attention due to their lower computational cost and favorable recovery performance. One class of methods parameterizes the unknown tensor through its core tensor and factor matrices, and then applies alternating minimization or gradient descent to the resulting factorized problem [28, 29, 30]. Another class of methods directly enforces the low-rank constraint through iterative hard thresholding (IHT) [15], which has been extended to the tensor setting as tensor iterative hard thresholding (TIHT) [21]. However, applying the truncated higher-order singular value decomposition (HOSVD) to a fulldimensional tensor requires computing the singular value decomposition of each large mode-i matricization, whose column dimension grows rapidly with the tensor order and can become prohibitively costly for high-dimensional problems. Riemannian gradient descent (RGD) overcomes this bottleneck by projecting the Euclidean gradient onto the tangent space of the fixed-multilinear-rank manifold at the current iterate, so that the retraction step only needs to process a compact core tensor rather than the full ambient tensor. Kressner et al. first developed this Riemannian optimization framework for tensor completion, introducing the tangent-space parametrization and HOSVD-based retraction on which subsequent RGD methods rely [31]. Cai et al. extended RGD to general linear measurement operators and proved that, once initialized by one step of IHT, the resulting iteration converges linearly to the underlying tensor under a tensor restricted isometry condition, with a sampling complexity optimal in the ambient dimension n and every retraction confined to core tensors of size at most 2r1 × · · · × 2rd [32]. More recent work has further extended RGD along several directions, including its implicit regularization behavior [33] and variants based on preconditioned or adapted Riemannian metrics [34, 35, 36]. The theoretical recovery guarantees of the aforementioned methods rely on the TRIP. Existing TRIP results, however, are predominantly established for vectorized measurement operators. For a tensor in Q Rn1 ×···×nd , such an operator is represented by a matrix of size M × di=1 ni . The storage of this matrix can exceed that of the original tensor and therefore limits the practical use of vectorization-based measurements in large scale settings. To address this difficulty, Haselby et al. introduced modewise measurement operators based on tensor reshaping and mode products [18]. Compared with vectorized measurements, modewise operators require fewer random parameters, reduce the storage cost of the sensing map, and are compatible with structured and parallel implementations. Building upon the modewise measurement framework in [18], subsequent work extended this low-memory approach to Tucker approximation from one-pass streamed measurements [37]. The RGD and modewise measurement frameworks address two complementary limitations of large-scale tensor recovery. RGD reduces the computational cost of the optimization step by exploiting the geometry of the fixed-rank tensor manifold, whereas modewise operators reduce the storage cost of the measurement process. Despite its computational benefits, standard RGD suffers from a limitation, it applies a uniform global scaling to all orthogonal components of the Riemannian gradient. This approach neglects the heterogeneous scale disparities that often emerge between the core tensor and the factor matrices during optimization. Treating these distinct components with a single step size can lead to severe update imbalances, such as overshooting in the factors or stagnating in the core which degrades convergence speed. To mitigate such ill-conditioning, recent studies have explored Riemannian preconditioning and adapted metrics for low-rank matrix and tensor optimization [36, 38, 34, 35]. Building upon the rationale of these geometric refinements, our motivation is to design a computationally lightweight yet effective scaling strategy that dynamically corrects these imbalances. To this end, we propose the adaptive block-weighted method, which adaptively rescales the core and factor gradient components of the Riemannian gradient according to their individual magnitudes, thereby harmonizing their updates and accelerating the overall recovery process. Contributions In this study, we develop an adaptive block-weighted modewise Riemannian gradient descent method for low-Tucker-rank tensor recovery. The proposed approach combines the storage advantages of structured modewise measurements with the geometric efficiency of Riemannian optimization on the fixed-multilinear-rank tensor manifold. The main contributions are summarized as follows.

• We incorporate one-stage and two-stage modewise measurement operators into the RGD framework. The resulting recovery model avoids the storage of a single dense vectorized measurement matrix and remains compatible with modewise operators satisfying the tensor restricted isometry property. • We propose a normalized adaptive block-weighted tangent operator that rescales the core and factor components of the Riemannian gradient according to their regularized Frobenius norms. The normalization maintains a unit average of the block weights and recovers the standard RGD direction when the Riemannian gradient components have equal Frobenius norms. Moreover, the weighting does not increase the multilinear-rank bound of the search direction or the size of the reduced core tensor used in the retraction. • We establish a local linear convergence guarantee that quantifies the effect of the adaptive weighting. Under a TRIP(δ2r , 2r) condition and a sufficiently accurate initialization, the iterates remain in a prescribed neighborhood of the underlying tensor and converge to it at a linear rate for every admissible weight deviation. Paper Outline. The remainder of this paper is organized as follows. In Section 2, we introduce the tensor notation, Tucker decomposition, fixed-rank tensor manifold, Riemannian gradient descent framework, and modewise measurement operators. In Section 3, we present the proposed adaptive block-weighted modewise Riemannian gradient descent method together with the normalized weighting strategy and its efficient implementation. In Section 4, we establish the modewise TRIP conditions, local linear convergence guarantees, and sampling complexity of the proposed method. Section 5 presents numerical experiments under Gaussian and SORS measurements to demonstrate the recovery performance and computational efficiency of our method. Finally, Section 6 concludes the paper and discusses directions for future work.

Preliminaries and Problem Formulation

Throughout this paper, calligraphic letters, such as X , are used to denote tensors, while capital Roman letters, such as V and A, denote matrices. Bold lowercase letters, such as x, denote vectors. Linear operators on tensors are denoted by bold script letters, such as L . Sets and spaces are denoted by calligraphic or blackboard bold letters, depending on the context. For a tensor X ∈ Rn1 ×···×nd , its (j1 , . . . , jd )-th entry is denoted by xj1 ···jd or [X ]j1 ,...,jd . For any positive integer p, we write [p] := {1, . . . , p}; in particular, [d] denotes the set of tensor-mode indices. Integer multi-indices are denoted by bold lowercase letters, such as a = (a1 , . . . , ad ). For two multi-indices a, b ∈ Nd , the relation a ⪯ b means that ai ≤ bi for every i ∈ [d]. For a matrix X, its singular values are arranged in nonincreasing order, and σj (X) denotes its j-th largest singular value. The tensor inner product and the associated Frobenius norm are denoted by ⟨·, ·⟩F and ∥ · ∥F , respectively. The identity operator is denoted by I , the adjoint of a linear operator L is denoted by L ∗ , and its operator norm is denoted by ∥L ∥. In the remainder of this section, we introduce the basic tensor operations, the geometry of fixed-multilinearrank tensors, and the measurement operators used throughout the paper.

2.1

Tensor notation

Pn1 Pnd For tensors X , Y ∈ Rn1 ×···×nd , the inner product is ⟨X , Y⟩F := i1 =1 · · · id =1 xi1 ···id yi1 ···id and p Frobenius norm is ∥X ∥F := ⟨X , X ⟩F . n1 ×···×nd is defined as a linear mapping that unfolds The mode-i matricization of a tensor X ∈ R Q n ×N X into a matrix X(i) ∈ R i i , where Ni := j̸=i nj . Specifically, this operation reshapes X such that each mode-i fiber Xj1 ,...,ji−1 ,:,ji+1 ,...,jd ∈ Rni constitutes a column of X(i) . Furthermore, we define Qd

the vectorization operator vec : Rn1 ×···×nd → R j=1 nj as the canonical linear bijection that maps the tensor to a column vector, with its inverse denoted by tenp1 ,...,ps := vec−1 for a target space Rp1 ×···×ps . 4

For a tensor X ∈ Rn1 ×···×nd and a matrix U ∈ Rpi ×ni , the mode-i product, denoted by X ×i U , yields a tensor in Rn1 ×···×ni−1 ×pi ×ni+1 ×···×nd . Its entries are explicitly defined by X ×i U j1 ,...,j

i−1 ,ℓ,ji+1 ,...,jd

ni X

xj1 ,...,ji ,...,jd Uℓ,ji .

ji =1

for all (j1 , . . . , ji−1 , ℓ, ji+1 , . . . , jd ) ∈ [n1 ] × · · · × [ni−1 ] × [pi ] × [ni+1 ] × · · · × [nd ]. When applied to a rank-one tensor Y = ⃝dj=1 v (j) composed of vectors v (j) ∈ Rnj , the mode-i product acts exclusively on the i-th factor vector, yielding (j) Y ×i U = ⃝i−1 ◦ U v (i) ◦ ⃝dj=i+1 v (j) . j=1 v

2.2

Tucker decomposition and fixed-rank tensor manifold

Let r = (r1 , . . . , rd ) ∈ Nd with 1 ≤ ri ≤ ni for all i ∈ [d]. The multilinear rank of a tensor X ∈ n ×···×n 1 d R is defined by mulrank(X ) := rank(X(1) ), . . . , rank(X(d) ) . We say that X has multilinear rank at most rN if mulrank(X ) ⪯ r. Equivalently, there exist subspaces Ui ⊆ Rni and dim(Ui ) = ri such that X ∈ di=1 Ui . (i) (i) For each i ∈ [d], let {v1 , . . . , vri } be an orthonormal basis of Ui , and define the corresponding (i) (i) factor matrix V (i) := v1 , . . . , vri ∈ Rni ×ri . Then there exists a core tensor B ∈ Rr1 ×···×rd such that X = B ×1 V

(1)

×2 · · · ×d V

(d)

r1 X k1 =1

···

rd X

(i)

B(k1 , . . . , kd ) ⃝di=1 vki .

(4)

kd =1

This representation is called an orthogonal Tucker decomposition of X . The factor matrices satisfy T V (i) V (i) = Iri , for i ∈ [d]. Moreover, mulrank(X ) = r if and only if the core tensor B has full multilinear rank r. The collection of all such tensors forms a smooth embedded submanifold in the ambient Euclidean tensor space, denoted by Mr := X ∈ Rn1 ×···×nd : mulrank(X ) = r , P Q whose dimension is given by dim(Mr ) = di=1 ri + di=1 (ni ri − ri2 ) [31, 32]. By differentiating the factors of the orthogonal Tucker decomposition X = B ×i∈[d] V (i) ∈ Mr and imposing the gauge conditions (V (i) )T V̇ (i) = 0 to eliminate the infinitesimal nonuniqueness of the representation, the tangent space at X is characterized as ( ) d X TX Mr = Ḃ ×i∈[d] V (i) + B ×j∈[d]\{k} V (j) ×k V̇ (k) Ḃ ∈ Rr1 ×···×rd , (V (k) )T V̇ (k) = 0 . k=1

Under these gauge conditions, the tangent space naturally admits a mutually orthogonal direct sum decomposition: (0) (1) (d) TX Mr = SX ⊕ SX ⊕ · · · ⊕ SX (5) (0)

with respect to the Frobenius inner product. Specifically, this consists of one core-variation block SX := (k) {Ḃ ×i∈[d] V (i) | Ḃ ∈ Rr1 ×···×rd } and d factor-variation blocks SX := {B ×j∈[d]\{k} V (j) ×k V̇ (k) | V̇ (k) ∈ Rnk ×rk , (V (k) )T V̇ (k) = 0} for each k ∈ [d]. (k) (k) For k = 0, . . . , d, let ΠX : Rn1 ×···×nd −→ SX denote the orthogonal projector onto the k(k) th tangent block SX . Then the orthogonal projector from the ambient tensor space onto TX Mr can Pd (k) be written as P TX Mr = k=0 ΠX . Importantly, any tangent vector ξ ∈ TX Mr inherits a strictly bounded structural complexity from this subspace configuration, yielding a multilinear rank bounded componentwise by 2r, i.e., mulrank(ξ) ⪯ 2r [31]. 5

2.3

Riemannian gradient descent algorithm

To minimize the least-squares objective f (X ) = 12 ∥L (X ) − y∥22 over the fixed-rank manifold Mr , the standard RGD algorithm iteratively updates the tensor via gradient projection and manifold retraction [32]. At the l-th iterate Xl ∈ Mr , the Euclidean gradient is given by Gl := ∇f (Xl ) = L ∗ L (Xl ) − y . Under the Frobenius metric, the Riemannian gradient is uniquely determined by projecting the Euclidean gradient onto the current tangent space Sl := TXl Mr , yielding P Sl Gl . Moving along the negative Riemannian gradient direction, the exact line-search step size αl is computed as ∥P Sl Gl ∥2F αl = argmin f Xl − αP Sl Gl = , ∥L (P Sl Gl )∥22 α≥0 provided that P Sl Gl ̸= 0. The equality follows directly from the orthogonality of the projector P Sl , which yields ⟨Gl , P Sl Gl ⟩F = ∥P Sl Gl ∥2F . Because the intermediate tensor Xl −αl P Sl Gl resides in the tangent space Sl and typically possesses a multilinear rank up to 2r, a retraction operator is required to pull it back onto the prescribed lowrank manifold Mr . A computationally efficient and widely adopted retraction is the truncated HOSVD operator H r , which retains the leading ri left singular vectors of each mode-i matricization [39]. The resulting truncation procedure is summarized in Algorithm 1. Algorithm 1 Truncated HOSVD H r , [40, 11] Require: Tensor Y ∈ Rn1 ×···×nd and target rank r = (r1 , . . . , rd ). 1: for i = 1, . . . , d do 2: Compute V (i) ∈ Rni ×ri as the ri dominant left singular vectors of Y(i) . 3: end for 4: B ← Y ×i∈[d] (V (i) )T . 5: return H r (Y) := B ×i∈[d] V (i) . Unlike the matrix case, computing best low-rank tensor approximations is generally NP-hard [41]. A practical alternative is the truncated HOSVD, which serves as a quasi-projection onto the multilinearrank-r manifold Mr [40, 11]. Specifically, letting P Mr denote the exact projection onto Mr , we have √ ∥Y − H r (Y)∥F ≤ d ∥Y − P Mr (Y)∥F , Y ∈ Rn1 ×···×nd . (6) Consequently, initialized by X0 = H r (L ∗ y), the standard RGD iteration is formulated as Xl+1 = H r (Xl − αl P Sl Gl ). This formulation allows the retraction to be efficiently evaluated on a compact core tensor of size at most 2r1 × · · · × 2rd , completely avoiding the HOSVD of the full-dimensional ambient tensor.

2.4

Restricted isometry properties

We first recall the restricted isometry property (RIP) on a prescribed set, which will be used to characterize the component matrices of the modewise measurement operators introduced later. The definition applies to an arbitrary subset of a normed vector space and is therefore not restricted to a particular low-dimensional model. Definition 2.1 (RIP(ε, S) property). Let S be a subset of a normed vector space and let A be a linear map. For 0 < ε < 1, we say that A satisfies the RIP(ε, S) property if (1 − ε)∥z∥2 ≤ ∥A (z)∥2 ≤ (1 + ε)∥z∥2 ,

z ∈ S.

(7)

Thus, the RIP requires the linear map to approximately preserve the norm of every element in the prescribed set S. For tensor recovery, the corresponding norm-preservation property is TRIP. In this case, the set consists of tensors whose multilinear rank is bounded by a prescribed rank tuple. 6

Definition 2.2 (TRIP(δ, r) property, [21]). Let r = (r1 , . . . , rd ) ∈ Nd and let 0 < δ < 1. A linear map A is said to satisfy the TRIP(δ, r) property if (1 − δ)∥X ∥2F ≤ ∥A (X )∥2F ≤ (1 + δ)∥X ∥2F for every tensor X ∈ Rn1 ×···×nd satisfying mulrank(X ) ⪯ r. The TRIP therefore requires the measurement operator to act as an approximate isometry on the set of tensors with multilinear rank bounded by r. In the following subsection, we introduce the modewise measurement operators and specify restricted isometry conditions on their component matrices that ensure the corresponding TRIP.

2.5

Modewise measurement

Following the multi-stage modewise measurement framework (3) introduced in [18], a modewise measurement operator is a structured linear map constructed by alternating tensor reshaping and modewise matrix multiplications. Before describing the corresponding measurement construction, we introduce two classes of structured sets. Suppose that X ∈ Rn×···×n , let q ≥ 2 be an integer dividing d, and set d′ := d/q. q

Definition 2.3 (The set S1,2 , [18]). Consider a set of vectors in Rn , let   q O  q S1 := v (j) : v (j) ∈ Sn−1 , j ∈ [q] ⊆ Sn −1 ,   j=1

and define S2 :=

x+y ∥x+y∥2

q x, y ∈ S1 , s.t. ⟨x, y⟩ = 0 , then set S1,2 := S1 ∪ S2 ⊆ Sn −1 .

The HOSVD requires the factor vectors within each mode to be mutually orthogonal. After modewise compression, however, the transformed factor vectors are generally only approximately orthogonal. The following set is therefore used to characterize the intermediate tensors produced by the first-stage compression. Definition 2.4 (Nearly orthogonal tensors, [18]). Let R ≥ 1, µ ≥ 0, θ ≥ 0, and let re = (e r1 , . . . , red′ ) ∈ ′ d m ×···×m ′ 1 d N . The set NR,µ,θ,er consists of all tensors Y ∈ R that admit the form Y=

re1 X k1 =1

···

red′ X

′

(i)

C(k1 , . . . , kd′ ) ⃝di=1 vki ,

kd′ =1 (i)

where C ∈ Rre1 ×···×erd′ is the core tensor and vki ∈ Rmi , ki ∈ [e ri ], i ∈ [d′ ], are the factor vectors. The factor vectors are not required to be mutually orthogonal or normalized. Instead, the core tensor and factor vectors satisfy the following conditions: (a)

(i) 2 ≤ R for every i ∈ [d′ ] and ki ∈ [e ri ]; 2

vki

(i)

(b) ⟨vki , vk′ ⟩ ≤ µ for every i ∈ [d′ ] and any ki ̸= ki′ ; i

(c) ∥C∥F = 1; (d) The mode-i subtensors of C are mutually orthogonal for every i ∈ [d′ ]. More precisely, for any distinct indices p ̸= p′ in [e ri ], we have re1 X k1 =1

rei−1

···

rei+1

ki−1 =1 ki+1 =1

···

red′ X

C(k1 , . . . , ki−1 , p, ki+1 , . . . , kd′ )C(k1 , . . . , ki−1 , p′ , ki+1 , . . . , kd′ ) = 0.

kd′ =1

(e) ∥Y∥F ≥ θ. The first stage groups every q consecutive modes of the original tensor and applies a separate measurement matrix along each reshaped mode. The resulting tensor is then vectorized and, when further dimension reduction is required, compressed by a secondary measurement matrix. Specifically, we define L 1 := vec ◦A ◦ R and L 2 := A2nd ◦ L 1 , where R is the tensor reshaping operator, A is the first stage modewise compression map, and A2nd is the second stage measurement matrix. We refer to L 1 and L 2 as the one-stage and two-stage modewise measurement operators, respectively. The complete construction can be summarized as R

vec

2nd X −−→ R(X ) −−→ A (R(X )) −−−→ L 1 (X ) −−− −→ L 2 (X ).

The first stage consists of tensor reshaping and modewise compression, whereas the optional second stage applies an additional vectorized compression to the first-stage output. The details of these operations are given below. 2.5.1

Step 1: tensor reshaping

For simplicity, suppose that ni = n and ri = r for all i ∈ [d]. Consider the reshaping operator d O

′

R −→

d O

s=1

i=1

which groups every q consecutive modes of a d-mode tensor into a single mode. Thus, R reduces the tensor order from d to d′ = d/q. R is the unique linear operator whose action on a rank-one tensor is given by   qs O ′ ′ e(s) , x(ℓ)  =: ⃝ds=1 x (8) R ⃝di=1 x(i) := ⃝ds=1  ℓ=1+q(s−1) q

e(s) ∈ Rn . where ◦ and ⊗ denote the outer product and the Kronecker product, respectively, andNv ′ q Suppose that X admits the Tucker decomposition (4). Then its reshaping Xe := R(X ) ∈ ds=1 Rn can be written as rq rq X X (s) e 1 , . . . , jd′ ) ⃝d′ v e B(j (9) X = ··· s=1 ejs , j1 =1

jd′ =1

where mulrank Xe ⪯ re := (rq , . . . , rq ), Be is obtained by reshaping the original core tensor B, and

(s)

ejs is the Kronecker product of the corresponding q Tucker factor vectors of X . Since Kronecker each v (s)

products of orthonormal factor vectors remain orthonormal, the vectors {e vjs }rjs =1 form an orthonormal ′ family for each s ∈ [d ]. 2.5.2

Step 2: first-stage modewise compression

The second step performs the first-stage compression by applying a separate measurement matrix along q each mode of the reshaped tensor. For every i ∈ [d′ ], let Ai ∈ Rmi ×n , mi ≪ nq . Define linear map q q A : Rn ×···×n −→ Rm1 ×···×md′ by A (Y) := Y ×1 A1 ×2 · · · ×d′ Ad′ .

(10)

Thus, A acts modewise on d′ -mode tensors. The corresponding first-stage tensor-valued measurement map on the original tensor space is F (X ) := A (R(X )). Let Xe := R(X ) be the reshaped tensor obtained in Step 1. By (9), using the action of a mode product on rank-one tensors, we obtain F (X ) = A (Xe) =

re1 X j1 =1

···

red′ X jd′ =1

(i) e 1 , . . . , jd′ ) ⃝d′ e B(j A v . i ji i=1

(11)

For compatibility with the vector-valued measurement model used in the recovery problem, we define the one-stage modewise measurement operator by vectorizing the first-stage output: L 1 (X ) := vec F (X ) = vec A (R(X )) , Q′ where L 1 (X ) ∈ RM1st , M1st := di=1 mi . When mi = m for all i ∈ [d′ ], the intermediate measure′ ment dimension is M1st = md . Let tenm1 ,...,md′ : RM1st −→ Rm1 ×···×md′ , denote the inverse vectorization operator. Since the adjoint of a mode-i multiplication by Ai is the mode-i multiplication by ATi , the adjoint of L 1 is L ∗1 (z) = R ∗ tenm1 ,...,md′ (z) ×1 AT1 ×2 · · · ×d′ ATd′ . The expansion (11) shows that the first-stage map acts on the reshaped factor vectors through the component matrices Ai . Accordingly, we impose RIP(ε, S1,2 ) on these matrices, where S1,2 is defined in Definition 2.3. The following proposition shows that these componentwise RIP conditions imply the TRIP for the first-stage modewise measurement operator. Proposition 2.5 (Theorem 3.3, [18]). Let r ≥ 2 and r := (r, . . . , r) ∈ Nd . Suppose that the reshaping operator R and the modewise map A are defined by (8) and (10), respectively. Assume that, for every q i ∈ [d′ ], the matrix Ai ∈ Rmi ×n satisfies the RIP(ε, S1,2 ) condition in (7). Define δ := 4d′ rd ε and suppose that δ < 1. Then the composite map A ◦ R satisfies the TRIP(δ, r) condition. More precisely, (1 − δ)∥X ∥2F ≤ ∥A (R(X ))∥22 ≤ (1 + δ)∥X ∥2F for every X ∈ Rn×···×n satisfying mulrank(X ) ⪯ r. Moreover, since vectorization preserves the Frobenius norm, the one-stage vector-valued operator L 1 = vec ◦A ◦ R satisfies the equivalent bound (1 − δ)∥X ∥2F ≤ ∥L 1 (X )∥22 ≤ (1 + δ)∥X ∥2F . 2.5.3

Step 3: secondary compression

The first stage output F (X ) = A (R(X )) ∈ Rm1 ×···×md′ retains a tensor structure. However, the n o (i) rei ej transformed factor vectors Ai v , i ∈ [d′ ] appearing in (11) are generally no longer mutually j=1

orthogonal. Consequently, (11) need not be an HOSVD of F (X ), and its multilinear rank cannot be inferred directly from that representation. This difficulty is addressed using the nearly orthogonal tensor class introduced in Definition 2.4. For the theoretical statements below, we assume that mi = m for all ′ i ∈ [d′ ], so that M1st = md . ′

Proposition 2.6 (Lemma 5.11, [18]). Let r = (r, . . . , r) ∈ Nd and re = (rq , . . . , rq ) ∈ Nd . Suppose that X ∈ Rn×···×n is a unit-norm tensor (i.e., ∥X ∥F = 1) satisfying mulrank(X ) ⪯ r. Assume that the first-stage modewise map A is defined by (10), and each matrix Ai satisfies the RIP(ε, S1,2 ) property for i ∈ [d′ ]. If δ := 12d′ rd ε < 1, then the first-stage output satisfies: A (R(X )) ∈ N1+ε,ε,1−δ/3,er . Proposition 2.6 demonstrates that the first-stage measurements of all unit-norm tensors with multilinear rank bounded by r reside in the set of nearly orthogonal tensors N1+ε,ε,1−δ/3,er . Consequently, this specific set defines the domain over which the secondary compression matrix must preserve Euclidean norms. To achieve a more compact final measurement vector, we introduce a secondary measurement matrix A2nd ∈ RM2nd ×M1st . The cascaded two-stage modewise measurement operator is then defined as: L 2 (X ) := A2nd L 1 (X ) = A2nd vec A (R(X )) ∈ RM2nd , (12) with adjoint given by L ∗2 (z) = L ∗1 (AT2nd z) for any z ∈ RM2nd . The following proposition gives a sufficient condition for the two-stage operator to satisfy the TRIP. 9

Proposition 2.7 (Theorem 3.8, [18]). Let r ≥ 2, r = (r, . . . , r) ∈ Nd , and re = (rq , . . . , rq ) ∈ ′ Nd . Suppose R, A , and L 2 are defined by (8), (10), and (12), respectively. Assume that each Ai ∈ q Rmi ×n satisfies the RIP(ε, S1,2 ) property for i ∈ [d′ ], and let δ = 12d′ rd ε < 1. If A2nd satisfies the RIP δ/3, vec(N1+ε,ε,1−δ/3,er ) property, then L 2 satisfies the TRIP(δ, r) condition. Specifically, (1 − δ)∥X ∥2F ≤ ∥L 2 (X )∥22 ≤ (1 + δ)∥X ∥2F holds for every X ∈ Rn×···×n satisfying mulrank(X ) ⪯ r. Proposition 2.7 shows that, under the prescribed RIP conditions, the secondary compression further reduces the measurement dimension while preserving the approximate isometry of the modewise measurement operator over tensors of bounded multilinear rank. In Section 4.1, we specialize these guarantees to multilinear rank bounded by 2r and present the corresponding measurement dimension bounds for sub-Gaussian and SORS matrices. In the remainder of this paper, L denotes either L 1 or L 2 , depending on whether the one-stage or two-stage measurement model is utilized.

Adaptive Block-Weighted Modewise RGD

To recover the underlying low-rank tensor from the compressed measurements, we formulate the task as a constrained least squares optimization problem over the fixed rank tensor manifold Mr : min f (X ) :=

X ∈Mr

1 ∥L (X ) − y∥22 . 2

(13)

In this section, we introduce the Adaptive Block-Weighted Modewise Riemannian Gradient Descent algorithm for solving the low-rank tensor recovery problem in (13).

3.1

Adaptive weighted tangent operator and algorithm

We construct the adaptive weighted tangent operator from the gradient at each iteration. Let Gl := L ∗ L (Xl ) − y be the Euclidean gradient of the objective function at the l-th iterate Xl , and let Sl := TXl Mr denote the tangent space of the manifold Mr at Xl . In our proposed algorithm, we retain the standard orthogonal tangent space decomposition and adaptively rescale its individual blocks. As introduced in Section 2.2, the tangent space Sl admits the mutually orthogonal direct-sum de(0) (1) (d) (k) (k) composition Sl = Sl ⊕ Sl ⊕ · · · ⊕ Sl , where Πl denotes the orthogonal projector onto Sl . P (k) Consequently, the orthogonal projector onto the full tangent space can be written as P Sl = dk=0 Πl . By applying these block projectors to the Euclidean gradient Gl , we define the gradient components as follows: (0) (k) (k) Dl = Πl Gl , Wl = Πl Gl , k = 1, . . . , d. (14) This yields the standard decomposition of the Riemannian gradient: P Sl Gl = Dl +

d X

(k)

Wl .

(15)

k=1

In standard RGD, all components in (15) contribute to the Riemannian gradient with the same unit (k) weight. However, the relative magnitudes of the core and factor gradient components, given by Πl Gl for k = 0, . . . , d, can vary considerably across iterations. Motivated by this observation, we introduce an adaptive weighting rule based on the Frobenius norms of these individual tangent-gradient components. Specifically, we define the block magnitudes as r sl,k =

(k)

Πl Gl

+ η2, 10

k = 0, . . . , d,

(16)

where η > 0 is a strictly positive regularization parameter. The regularization parameter η > 0 keeps all regularized block magnitudes strictly positive and ensures that, as the tangent-gradient components vanish near convergence, the normalized weights approach one, thereby recovering the standard RGD scaling. To convert these regularized block magnitudes into comparable scaling factors for the tangent components, we normalize them so that the average weight remains equal to one. This leads to the following definition. Definition 3.1 (Normalized adaptive weights). At the l-th iteration, the normalized adaptive weight corresponding to the k-th tangent block is defined as sl,k ωl,k = (d + 1) Pd

j=0 sl,j

k = 0, . . . , d.

(17)

The normalization in Definition 3.1 preserves the average scale of the tangent-gradient blocks. The weighting mechanism redistributes the relative contributions of the tangent blocks without introducing an additional iteration-dependent global scaling. The following elementary bound will be useful in Section 4. (k)

Lemma 3.2 (Boundedness of the adaptive weights). Suppose that max0≤k≤d ∥Πl Gl ∥F ≤ Γ on a set of iterates. Then, for all k = 0, . . . , d, η 0< p =: ωmin ≤ ωl,k ≤ d + 1 =: ωmax . 2 Γ + η2 Based on the adaptive weights defined above, we obtain an adaptive block-weighted tangent operator as follows. Definition 3.3 (Adaptive block-weighted tangent operator). Let {ωl,k }dk=0 be the normalized adaptive weights defined in Definition 3.1. The adaptive block-weighted tangent operator at Xl is defined as l Pω l =

d X

(k)

ωl,k Πl ,

(18)

k=0

The corresponding adaptive block-weighted tangent direction is defined by l Zl = P ω l Gl .

(19)

Remark 3.4. If all tangent-gradient components have the same Frobenius norm, then Definition 3.1 gives l ωl,0 = · · · = ωl,d = 1. Hence P ω l = P Sl , and Zl = P Sl Gl . Thus, the proposed weighting scheme recovers the standard RGD search direction as a special case. Using (14), the weighted tangent gradient admits the block representation Zl = ωl,0 Dl +

d X

(k)

ωl,k Wl .

k=1

The adaptive operator (18) also admits a variable-metric interpretation [38]. For ξ, ν ∈ Sl , define the iteration-dependent inner product glωl (ξ, ν) =

d X

D E (k) (k) −1 ωl,k Πl ξ, Πl ν .

(20)

k=0

Since ωl,k > 0 for every k, this defines a positive definite inner product on Sl . For any ξ ∈ Sl , using (19)–(20), we obtain glωl (Zl , ξ) =

d X k=0

−1 ωl,k

(k) (k) Πl Zl , Πl ξ

E F

d D X k=0

(k)

Πl Gl , Πl ξ

E F

= ⟨Gl , ξ⟩F = Df (Xl )[ξ],

Algorithm 2 Adaptive Block-Weighted Modewise RGD 1: Initialize: X0 = H r L ∗ (y) 2: for l = 0, 1, . . . do 3: Compute gradient: Gl = L ∗ L (Xl ) − y 4: Update adaptive weights: Compute ωl,k for k = 0, 1, . . . , d according to (17) l 5: Compute search direction: Compute Zl = P ω l Gl according to (18) and (19) 6: Compute step size: Compute αl according to (21) 7: Update tensor: Xl+1 = H r Xl − αl Zl 8: end for 9: Output: Xl+1 when the stopping criterion is met. where Df (Xl )[ξ] denotes the directional derivative of the objective function f at Xl along the tangent direction ξ. Therefore, Zl can be viewed as the gradient associated with the iteration-dependent metric (20). Since the weights depend on the current gradient, this metric is understood as a local variable metric defined at each iteration rather than a fixed metric prescribed globally on Mr . Having constructed the adaptive tangent direction Zl , we next determine the step size by minimizing the least-squares objective along the search direction −Zl . Since the objective is quadratic in α, the exact line search admits a closed-form solution. Specifically, d X

αl = argmin f (Xl − αZl ) = α≥0

⟨Gl , Zl ⟩F

(k)

ωl,k Πl Gl

= k=0 ∥L (Zl )∥22 ∥L (Zl )∥22

2 F

(21)

provided that L (Zl ) ̸= 0. In general, the numerator in (21) is not equal to ∥Zl ∥2F . By the orthogonality of the tangent blocks, 2 Pd (k) 2 . The equality ⟨Gl , Zl ⟩F = ∥Zl ∥2F holds only in the standard un∥Zl ∥2F = k=0 ωl,k Πl Gl F

l weighted case. This distinction is important because P ω l is not an orthogonal projection and will also play a role in the convergence analysis in Section 4. Building on the adaptive block-weighted tangent operator and the associated variable metric introduced above, we formulate the block-weighted modewise RGD method for solving (13). At each iteration, the residual is first mapped back to the ambient tensor space by L ∗ to obtain the Euclidean gradient, which is then decomposed into its core and factor tangent components. The normalized adaptive weights are updated from the Frobenius norms of these components, and the corresponding block-weighted tangent direction is constructed according to Definition 3.3. An exact line search is subsequently performed along this direction, followed by a truncated-HOSVD retraction onto Mr . The complete procedure is summarized in Algorithm 2.

3.2

Computation in Adaptive Block-Weighted RGD

3.2.1

Computation of the search direction

We need to compute the d + 1 tangent-gradient components introduced in (14). Because the direct sum decomposition of the tangent space Sl in (5) is mutually orthogonal, the orthogonal projection onto Sl can be decomposed into independent least-squares subproblems for the core tensor component and each factor matrix component. Let Y ∈ Rn1 ×···×nd be an arbitrary tensor, its orthogonal projection onto Sl is given by the solution to the following optimization problem: P Sl Y = argmin ∥Y − ξ∥2F . ξ∈Sl

(22)

An arbitrary element ξ ∈ Sl admits the representation (i)

ξ = Ḃ ×i∈[d] Vl

d X

(j)

Bl ×j∈[d]\{k} Vl

×k V̇ (k) ,

(23)

k=1 (k)

where Vl

(k) T , V̇ (k) ∈ Rnk ×rk satisfy Vl V̇ (k) = 0 for k = 1, . . . , d. (i)

The d + 1 terms in (23), consisting of the core-variation term Ḃ ×i∈[d] Vl

and the d factor-variation

(j) terms Bl ×j∈[d]\{k} Vl ×k V̇ (k) , k

= 1, . . . , d, are mutually orthogonal. Consequently, the global least-squares problem (22) decouples into d + 1 independent subproblems [32]. Core Component

The core variation Ḃl is obtained by solving Ḃl =

(i) 2

Y − Ḃ ×i∈[d] Vl

argmin Ḃ∈Rr1 ×···×rd

(j) T The closed-form solution is Ḃl = Y ×dj=1 Vl . Consequently, the projection onto the core tangent (0)

(j)

subspace is Πl Y = Ḃl ×dj=1 Vl . Factor Components the k-th mode) as

For k ∈ [d], we define the Kronecker product of the factor matrices (excluding (k)

V̄l

(d)

= Vl

(k+1)

⊗ · · · ⊗ Vl

(k−1)

⊗ Vl

(1)

⊗ · · · ⊗ Vl

The k-th factor variation is obtained by solving the constrained matrix least-squares problem (k)

V̇l

(k) T Y(k) − V̇ (k) (Bl )(k) V̄l

argmin (k) ∈Rnk ×rk V̇ (k) T (k) V̇ =0 Vl

. F

(k) (k) T Let P V (k)⊥ = Ink − Vl Vl denote the orthogonal projector onto the orthogonal complement of l (k) the column space of Vl . The closed-form solution is (k)

V̇l

(k)

= P V (k)⊥ Y(k) V̄l l

(Bl )†(k) ,

k = 1, . . . , d,

(24)

where Y(k) and (Bl )(k) denote the mode-k matricizations of Y and Bl , respectively, and (Bl )†(k) is the Moore–Penrose pseudoinverse of (Bl )(k) . The corresponding projection onto the k-th factor tangent (k)

(j)

block is Πl Y = Bl ×j∈[d]\{k} Vl

(k)

×k V̇l

. (0)

By setting Y = Gl , we obtain the core and factor tangent-gradient components Dl := Πl Gl = (i) (k) (k) (j) (k) Ḃl ×i∈[d] Vl and Wl := Πl Gl = Bl ×j∈[d]\{k} Vl ×k V̇l for k = 1, . . . , d. Consequently, the adaptive block-weighted tangent direction admits the representation (i)

Zl = ωl,0 Ḃl ×i∈[d] Vl

d X

(j)

ωl,k Bl ×j∈[d]\{k} Vl

(k)

×k V̇l

k=1

Remark 3.5 (Computation and scaling of the block magnitudes). The adaptive weights introduced in Definition 3.1 are determined by the Frobenius norms of the tangent-gradient blocks. These block norms can be efficiently evaluated without forming the corresponding full tangent tensors in the ambient space. (i) For the core component Dl = Ḃl ×i∈[d] Vl , exploiting the orthonormality of the factor matrices yields (k)

∥Dl ∥F = ∥Ḃl ∥F . For the factor components, the k-th component Wl 13

(j)

= Bl ×j∈[d]\{k} Vl

(k)

×k V̇l

has

(k)

the mode-k matricization (Wl )(k) = V̇l

(k) T

(Bl )(k) (V̄l

(k)

) . Because V̄l

has orthonormal columns,

(k) it follows that ∥Wl ∥F q

(k) = ∥V̇l (Bl )(k) ∥F . Thus, the block magnitudes in (16) can be computed as q (k) sl,0 = ∥Ḃl ∥2F + η 2 and sl,k = ∥V̇l (Bl )(k) ∥2F + η 2 , k = 1, . . . , d. This formulation avoids con(k) structing the full ambient tangent tensors when computing the adaptive weights. Moreover, Ḃl and V̇l

belong to different blocks of the Tucker parametrization, so their parameter-space norms are not directly (k) comparable. The quantities ∥Dl ∥F and ∥Wl ∥F instead measure the corresponding perturbations in the common ambient Frobenius norm. Since the adaptive weighting is applied after the orthogonal tangent-space projection, the least(k) squares subproblems for computing Ḃl and V̇l are unchanged. The additional work consists only of evaluating the d + 1 block magnitudes and normalizing the corresponding scalar weights. 3.2.2

Computation of the Retraction

We now detail the computation of the retraction step in Algorithm 2. Let Yl = Xl − αl Zl denote the intermediate tensor in the tangent space Sl . To obtain the next iterate, we retract Yl from the tangent space Sl back to the smooth multilinear-rank manifold Mr via the truncated HOSVD operator, i.e., Xl+1 = H r (Yl ). A direct implementation would require applying the truncated HOSVD to the full ambient tensor Yl . Using the tangent-space representation, the same retraction can instead be computed from a reduced core tensor [32]. Since both Xl and the search direction Zl reside in Sl , the updated tensor Yl belongs to Sl . Consequently, its multilinear rank is bounded by 2r (i.e., mulrank(Yl ) ⪯ 2r). Specifically, expanding the terms yields ! d X (k) (j) (i) ωl,k Bl ×j∈[d]\{k} Vl ×k V̇l Yl = Xl − αl Zl = Xl − αl ωl,0 Ḃl ×i∈[d] Vl + k=1

(i)

= Bl − αl ωl,0 Ḃl ×i∈[d] Vl

− αl

d X

(j)

ωl,k Bl ×j∈[d]\{k} Vl

(k)

×k V̇l

h i := Cl ×i∈[d] Vl(i) V̇l(i) .

k=1

Thus, Yl admits an augmented Tucker representation parameterized by a sparse block core tensor Cl ∈ R2r1 ×···×2rd . The principal block of this core tensor is given by [Cl ]1:r1 ,...,1:rd = Bl − αl ωl,0 Ḃl . Furthermore, for each k = 1, . . . , d, the block corresponding to the variation of the k-th factor matrix is [Cl ]1:r1 ,...,rk +1:2rk ,...,1:rd = −αl ωl,k Bl . The weights ωl,k scale only the nonzero blocks of Cl and do not introduce additional factor directions. Hence the multilinear-rank bound is unchanged. To compute H r (Yl ) efficiently, we first compute the thin QR factorizations of the augmented factor matrices: h i (i) (i) (i) (i) = Ql Rl , i = 1, . . . , d. Vl V̇l Using these orthogonal factors, Yl can be equivalently expressed as (i) (i) (i) (i) (i) Yl = Cl ×i∈[d] Ql Rl = Cl ×i∈[d] Rl ×i∈[d] Ql := Cel ×i∈[d] Ql . (i)

(i)

where Cel = Cl ×i∈[d] Rl ∈ R2r1 ×···×2rd . Since the matrices Ql possess orthonormal columns for all i = 1, . . . , d, the dominant mode subspaces of Yl exactly correspond to the leading singular vectors of the respective matricizations of the reduced core tensor Cel . Specifically, for each mode i = 1, . . . , d, let e (i) denote the matrix of left singular vectors of (Cel )(i) . The updated factor matrices are obtained via U l h i (i) (i) e (i) Vl+1 = Ql U , i = 1, . . . , d. l :,1:ri

The updated core tensor can be computed entirely in the reduced coordinates as h T i (i) e e Bl+1 = Cl ×i∈[d] Ul . :,1:ri

Consequently, the final retracted tensor is formulated as (i)

Xl+1 = Bl+1 ×i∈[d] Vl+1 = H r (Yl ). Thus, the truncated HOSVD need only be applied to the reduced tensor Cel ∈ R2r1 ×···×2rd , together with d thin QR factorizations of the augmented factor matrices. The adaptive weighting leaves the reduced-core retraction unchanged and adds only the evaluation of d + 1 block magnitudes and their scalar normalization. The measurement operator and its adjoint are applied through the modewise implementations of L and L ∗ . Hence the search direction retains multilinear rank at most 2r, and the truncated HOSVD is performed on a tensor of size 2r1 × · · · × 2rd .

Convergence Analysis

In this section, we establish the local linear convergence guarantee of Algorithm 2. In particular, we prove that Algorithm 2 converges linearly to the underlying unknown tensor T , provided that the modewise measurement operator satisfies the TRIP. Furthermore, we establish the sampling complexity of the proposed approach.

4.1

Modewise Tensor Restricted Isometry Property

For the explicit modewise TRIP and sampling-complexity bounds, we restrict to the balanced setting ni = n, ri = r for i ∈ [d]. We first specify the order of TRIP required for the subsequent theoretical Pd (k) l introduced in Secanalysis. In the adaptive block-weighted tangent operator P ω k=0 ωl,k Πl l = (k) tion 3, Πl is the orthogonal projector onto the k-th component of Sl . Since scalar multiplication by (k) (k) ωl,k does not alter the underlying subspace, then ωl,k Πl Gl ∈ Sl for k = 0, . . . , d. Consequently, the ωl weighted projected gradient satisfies Zl = P l Gl ∈ Sl with mulrank(Zl ) ⪯ 2r, which implies that the intermediate tensor Yl = Xl − αl Zl has a multilinear rank of at most 2r. The proposed adaptive weighting scheme does not increase the multilinear rank order required by the TRIP. Next, we describe the modewise measurement operators that satisfy this TRIP assumption. We first consider the one-stage modewise measurement operator L 1 introduced in Section 2.5, with modewise q measurement matrices Ai ∈ Rm×n , i = 1, . . . , d′ . The following proposition gives the corresponding TRIP guarantee for tensors of multilinear rank at most 2r. Proposition 4.1 (One-stage modewise TRIP at rank 2r, [18]). Suppose r ≥ 1, q ≥ 2, and 0 < ζ < 1, and define d′ := d/q. If each measurement matrix Ai satisfies the RIP(ε, S1,2 ) property, and δ2r = 4d′ (2r)d ε < 1, then the first-stage measurement map L 1 respects the TRIP(δ2r , 2r) property, i.e., (1 − δ2r )∥X ∥2F ≤ ∥L 1 (X )∥22 ≤ (1 + δ2r )∥X ∥2F holds for all tensors X with multilinear rank at most 2r. Furthermore, this conclusion holds with probability at least 1−ζ provided that either of the following measurement ensembles is utilized: (i) If the entries of each Ai are properly normalized i.i.d.sub-Gaussian random variables, it is sufficient that 2 nd log q d2 d −2 2d m ≥ Cδ2r (2r) max , 2 log . q q qζ (ii) If each Ai is a SORS matrix, it is sufficient that −2 m ≥ C1 δ2r (2r)2d

nd2 log q Ψ1 , q

where the logarithmic factor Ψ1 is defined as 2 2enq d 2d 2d −2 2 2d nd log q Ψ1 := log log log C2 δ2r (2r) log . qζ qζ q qζ 15

(25)

We next consider the two-stage modewise measurement operator L 2 introduced in Section 2.5, d′ with second-stage matrix A2nd ∈ RM2nd ×m . The following proposition gives the corresponding TRIP guarantee for tensors of multilinear rank at most 2r. Proposition 4.2 (Two-stage modewise TRIP at rank 2r, [18]). Suppose r ≥ 1, q ≥ 2, and 0 < ζ < 1, and define d′ := d/q. Assume that each first-stage measurement matrix Ai satisfies the RIP(ε, S1,2 ) property, and let δ2r = 12d′ (2r)d ε < 1. If the second-stage measurement matrix A2nd satisfies the associated RIP condition with reshaped rank ((2r)q , . . . , (2r)q ), then the two-stage measurement map L 2 respects the TRIP(δ2r , 2r) property, i.e., (1 − δ2r )∥X ∥2F ≤ ∥L 2 (X )∥22 ≤ (1 + δ2r )∥X ∥2F holds for all tensors X with multilinear rank at most 2r. Furthermore, this conclusion holds with probability at least 1−ζ provided that either of the following measurement ensembles is utilized: (i) If the measurement matrices consist of properly normalized i.i.d. sub-Gaussian random variables, it is sufficient that 2 nd log q d2 2d −2 2d m ≥ Cδ2r (2r) max , 2 log , q q qζ and −2 max M2nd ≥ Cδ2r

(

(2r)d q + dm(2r)q q

log

d +1 q

d2 m(2r)q δ dm(2r)q 2r + log 1 + δ2r (2r)d + , log 2 q q

) 2 . ζ

(ii) If the measurement matrices are SORS matrices, the first-stage requirement on m is identical to (25). For the second stage, it is sufficient that " (2r)d q + dm(2r)q d −2 M2nd ≥ Cδ2r log +1 q q # d2 m(2r)q δ dm(2r)q 2r log 1 + δ2r (2r)d + Ψ2 , (26) + q q2 where the logarithmic factor Ψ2 is defined as " c1 4 (2r)d q + dm(2r)q d 2 Ψ2 := log log +1 2 log ζ q q δ2r d2 m(2r)q δ dm(2r)q 2r + log 1 + δ2r (2r)d + 2 q q 4 4em × log log . ζ ζ

4.2

Main Theorems

We now establish the local linear convergence of the proposed adaptive block-weighted RGD iteration. To quantify the effect of the adaptive weighting, we first measure the deviation of the adaptive weights from the unit weights used in the standard RGD direction. Specifically, we define the maximal weight 16

deviation at the l-th iteration as ρl := max0≤k≤d |ωl,k − 1|. To bound this deviation, let Γ p ≥ 0 be an upper bound on the Frobenius norms of the tangent-gradient components, and define SΓ := Γ2 + η 2 . Then SΓ provides an upper bound on the corresponding regularized block magnitudes, and ρη (Γ) := d(SΓ −η) SΓ +dη gives a uniform upper bound on the deviation of the normalized adaptive weights from one. The weight deviation can further be controlled in terms of the current reconstruction error. As shown in Appendix 6, the TRIP(δ2r , 2r) condition yields an upper bound on the tangent-gradient components. Hence, whenever ∥Xl − T ∥F ≤ R, we have ρl ≤ ρη ((1 + δ2r )R). The auxiliary geometric estimates, the TRIP-based bounds, and the complete convergence proof are provided in Appendix 6. Theorem 4.3 (Local linear convergence of adaptive block-weighted modewise RGD). Let T ∈ Mr and y = L (T ). Assume that the measurement operator L satisfies TRIP(δ2r , 2r) with δ2r ∈ (0, 1). Fix a local radius R > 0 and an admissible weight-deviation level ρ⋆ ∈ (0, 1). Suppose that the initial iterate X0 ∈ Mr satisfies ∥X0 − T ∥F < R, and choose the regularization parameter η > 0 such that ρη ((1 + δ2r )R) ≤ ρ⋆ . Define √ 2( d + 1) (2d − 1)(1 + δ2r ρ⋆ ) γR := δ2r + ρ⋆ + R . (1 − δ2r )(1 − ρ⋆ ) min1≤i≤d σri (T(i) ) If γR < 1, then the iterates generated by Algorithm 2 satisfy ∥Xl − T ∥F ≤ γRl ∥X0 − T ∥F < R,

l ≥ 0.

Moreover, the adaptive weights satisfy ρl ≤ ρ⋆ for all l ≥ 0. Corollary 4.4 (Contraction under the standard RGD initialization). Suppose that the initialization satisfies √ (27) ∥X0 − T ∥F ≤ ( d + 1)δ2r ∥T ∥F . √ Fix ρ⋆ ∈ (0, 1) and choose η such that ρη (1 + δ2r )( d + 1)δ2r ∥T ∥F ≤ ρ⋆ . Then Theorem 4.3 √ applies with R = ( d + 1)δ2r ∥T ∥F . In particular, the contraction factor becomes " √ 2( d + 1) γw := δ2r + ρ⋆ (1 − δ2r )(1 − ρ⋆ ) # (28) √ ∥T ∥ F + (2d − 1)(1 + δ2r ρ⋆ )( d + 1)δ2r . min1≤i≤d σri (T(i) ) If γw < 1, then ∥Xl − T ∥F ≤ γwl ∥X0 − T ∥F . The following results further characterize the dependence of the local convergence bound on the weight deviation and the regularization parameter, as well as the asymptotic behavior of the adaptive weights. Remark 4.5 (Relation to the unweighted contraction factor). Setting ρ⋆ = 0 in (28) yields √ √ 2( d + 1)δ2r ∥T ∥F d γu := γw ρ⋆ =0 = 1 + (2 − 1)( d + 1) , 1 − δ2r min1≤i≤d σri (T(i) ) which agrees with the corresponding unweighted bound. Moreover, γw is increasing with respect to √ √ 2( d+1) N (ρ⋆ ) ρ⋆ on [0, 1). Indeed, writing γw = 1−δ2r 1−ρ⋆ , where N (ρ⋆ ) = δ2r + ρ⋆ + (2d − 1)( d + ∥T ∥F 1)δ2r min1≤i≤d σr (T i

∂ ∂ρ⋆

(i) )

(1 + δ2r ρ⋆ ), we have

N (ρ⋆ ) 1 − ρ⋆

h √ ∥T ∥F (1 + δ2r ) 1 + (2d − 1)( d + 1)δ2r min1≤i≤d σr (T i

(1 − ρ⋆ )2

Hence γw ≥ γu , with equality at ρ⋆ = 0. 17

i (i) )

> 0.

Remark 4.6 (Explicit choice of the regularization parameter). p The condition ρη ((1 +pδ2r )R) ≤ ρ⋆ in Theorem 4.3 can be ensured by an explicit choice of η. Since Γ2 + η 2 +η ≥ 2η and Γ2 + η 2 +dη ≥ (1 + d)η, (47) gives d Γ2 d Γ2 . ρη (Γ) = p p ≤ 2(1 + d)η 2 Γ2 + η 2 + η Γ2 + η 2 + dη Therefore, with Γ = (1 + δ2r )R, the sufficient condition η ≥ (1 + δ2r )R

(29)

d 2(1+d)ρ⋆ guarantees ρl ≤ ρ⋆

whenever ∥Xl − T ∥F ≤ R. Proposition 4.7 (Asymptotic weight deviation and contraction rate). Let the assumptions of Theorem 4.3 hold with γR < 1, so that Xl → T . Then the adaptive weights satisfy ρl ≤

d(1 + δ2r )2 ∥Xl − T ∥2F , 2(1 + d)η 2

(30)

that is, the weight deviation is second order in the reconstruction error, and the asymptotic contraction rate obeys √ ∥Xl+1 − T ∥F 2( d + 1)δ2r lim sup ≤ , (31) ∥Xl − T ∥F 1 − δ2r l→∞ which agrees with the corresponding asymptotic bound for the unweighted modewise RGD iteration. (k)

Proof. The block-gradient estimate (48) gives max0≤k≤d ∥Πl Gl ∥F ≤ (1 + δ2r )∥Xl − T ∥F =: Γl . Applying (29) with Γ = Γl yields (30). Since γR < 1 forces ∥Xl − T ∥F → 0, we also have ρl → 0. Substituting ρl → 0 and ∥Xl − T ∥F → 0 into the one-step recurrence (50) leaves the limiting factor √ 2( d + 1)δ2r /(1 − δ2r ), proving (31). Corollary 4.8 (Sampling complexity for linear convergence). Assume the balanced setting ni = n and ri = r for i ∈ [d]. Let q ≥ 2 divide d, d′ = d/q, ζ ∈ (0, 1) be the failure probability, and ρ⋆ ∈ (0, 1) be the target deviation level. Suppose δ⋆ ∈ (0, 1) satisfies √ √ ∥T ∥F d 2( d + 1) δ⋆ + ρ⋆ + (2 − 1)(1 + δ⋆ ρ⋆ )( d + 1)δ⋆ < (1 − δ⋆ )(1 − ρ⋆ ), (32) min1≤i≤d σri (T(i) ) √ and choose η > 0 such that ρη (1 + δ⋆ )( d + 1)δ⋆ ∥T ∥F ≤ ρ⋆ . Let y = L (T ) and X0 = H r (L ∗ y) be the initialization of Algorithm 2. Then the following hold: (i) For the one-stage sub-Gaussian modewise measurement map, if 2 nd log q d2 d −2 2d m ≥ Cδ⋆ (2r) max , 2 log , q q qζ then with probability at least 1 − ζ, Algorithm 2 converges linearly to T : ∥Xl − T ∥F ≤ γwl ∥X0 − T ∥F , ′

where γw < 1 is defined in (28). The total measurement dimension is M1st = md . (ii) For the two-stage sub-Gaussian modewise measurement map, if 2 nd log q d2 2d −2 2d m ≥ Cδ⋆ (2r) max , 2 log , q q qζ

and M2nd ≥ Cδ⋆−2 max

(

(2r)d q + dm(2r)q q

log

d +1 q

d2 m(2r)q δ dm(2r)q ⋆ log 1 + δ⋆ (2r)d + , log + 2 q q

) 2 , ζ

then with probability at least 1 − ζ, Algorithm 2 converges linearly to T with contraction factor γw < 1: ∥Xl − T ∥F ≤ γwl ∥X0 − T ∥F . Analogous convergence guarantees hold for SORS measurement ensembles by substituting the respective sampling conditions with (25) and (26). Proof. By Propositions 4.1 and 4.2, the stated sampling conditions ensure that L satisfies TRIP(δ⋆ , 2r) with probability at √least 1 − ζ. Under this TRIP condition, the truncated-HOSVD initialization satisfies ∥X0 − T ∥F ≤ ( d + 1)δ⋆ ∥T ∥F . The prescribed choice of η controls the adaptive weight deviation, while (32) ensures γw < 1. The conclusion therefore follows from Corollary 4.4. Remark 4.9 (Measurement dimension and rank dependence). For the one-stage construction, m is linear in n when the remaining parameters are fixed, but the total number of scalar outputs is M1st = md/q and is therefore superlinear in n whenever d/q > 1. In the two-stage construction, the final measurement dimension M2nd is linear in n up to the displayed rank and logarithmic factors.

Numerical Experiments

In this section, we evaluate the proposed normalized block-weighted modewise RGD method on synthetic low-multilinear-rank tensor recovery problems. The unknown fourth-order tensors are generated randomly in Tucker form. We consider both balanced and nonuniform tensor dimensions and multilinear ranks. The two-stage measurement operator is constructed from either Gaussian or SORS matrices. The intermediate modewise dimension is denoted by m, while p denotes the final target dimension. For comparison, we also test the unweighted modewise RGD method and the corresponding weighted and unweighted RGD methods with dense vectorized measurements. In the figures, “w” and “u” denote the weighted and unweighted variants, respectively, and “vec” denotes vectorized measurements. For each parameter setting, 20 independent trials are performed. A trial is regarded as successful if ∥Xl − T ∥F < 10−2 ∥T ∥F within 1000 iterations. The reported number of iterations is averaged over the 20 trials, with 1000 used as the iteration cap. Figures 1 and 2 report the results for Gaussian measurements. Four combinations of tensor dimensions and multilinear ranks are considered. In all cases, the successful recovery rate exhibits a clear transition as p increases. A larger intermediate dimension m generally shifts the transition to a smaller p and reduces the number of iterations. The weighted modewise method is comparable to, and often better than, its unweighted counterpart, especially near the recovery threshold and for smaller values of m. Figures 3 and 4 show the results for SORS measurements. In the balanced-rank cases, the modewise methods with m = 200 or 250 achieve nearly full recovery at target dimensions comparable to, and sometimes smaller than, those required by the vectorized methods. For the nonuniform rank (5, 6, 7, 8), increasing m from 150 to 200 or 250 substantially improves both the recovery rate and convergence. The weighted variant generally requires fewer iterations near the recovery threshold. To further evaluate the computational efficiency, Figure 5 compares five methods: the proposed weighted modewise RGD (weighted-mRGD), the unweighted modewise RGD (mRGD), the vectorized 19

Figure 1: Successful recovery rates for Gaussian measurements. Each point is computed from 20 independent trials, and recovery is declared successful when the relative error is below 10−2 within 1000 iterations.

Figure 2: Average iteration numbers for Gaussian measurements. Values close to 1000 correspond to parameter regimes in which many trials do not reach the prescribed accuracy.

Figure 3: Successful recovery rates for SORS measurements under balanced and unbalanced tensor dimensions and multilinear ranks.

Figure 4: Average iteration numbers for SORS measurements. Increasing the intermediate dimension improves convergence, particularly in the unbalanced-rank experiments.

Figure 5: Relative error versus CPU time for weighted-mRGD, mRGD, vec-RGD, TIHT, and vec-TIHT. The left panel uses Gaussian measurements with n = (20, 20, 20, 20) and r = (6, 6, 6, 6), while the right panel uses SORS measurements with n = (40, 40, 40, 40) and r = (6, 6, 6, 6). RGD (vec-RGD), the modewise TIHT (TIHT), and the vectorized TIHT (vec-TIHT), in terms of relative recovery error versus CPU time. For Gaussian measurements, weighted-mRGD exhibits the fastest convergence and reaches a relative error below 10−4 in approximately 7 seconds. Under SORS measurements, weighted-mRGD, mRGD, and vec-RGD converge, whereas TIHT and vec-TIHT remain nearly stagnant. Among the RGD methods, weighted-mRGD achieves the fastest error decay and the lowest final error. These results show that the normalized weighting improves the practical convergence efficiency of mRGD, particularly under SORS measurements. Figure 6 compares weighted-mRGD, mRGD, vec-RGD, TIHT, and vec-TIHT under Gaussian and SORS measurements. For Gaussian measurements, all methods achieve accurate recovery when p is sufficiently large, while weighted-mRGD generally requires fewer iterations near the recovery threshold. The difference is more evident for SORS measurements. Both weighted-mRGD and mRGD reach relative errors around 10−2 from p ≈ 6000, whereas vec-RGD requires a larger target dimension. In contrast, TIHT and vec-TIHT remain inaccurate and often reach the iteration cap over a large part of the tested range. These results show that combining modewise measurements with RGD improves the recovery efficiency, while the normalized weighting further accelerates the convergence of mRGD. Overall, the experiments confirm that the proposed method retains the recovery capability of vectorized RGD while using structured modewise measurements. The normalized adaptive weighting is most effective near the recovery threshold, where it generally improves the iteration count and, in several difficult settings, the empirical success rate.

Figure 6: Comparison of weighted and unweighted RGD and TIHT methods under Gaussian and SORS measurements. The top and bottom rows report the average iteration numbers and final relative errors, respectively. The left column corresponds to Gaussian measurements with n = (20, 20, 20, 20), r = (6, 6, 6, 6), and mint = 200, while the right column corresponds to SORS measurements with n = (40, 40, 40, 40), r = (6, 6, 6, 6), and mint = 250.

Conclusion and Future Work

In this paper, we combined modewise measurements with a normalized block-weighted Riemannian gradient framework for low-multilinear-rank tensor recovery. The proposed weighting adjusts the relative contributions of the tangent-gradient components while preserving the underlying multilinear-rank structure. Numerical results demonstrate improved convergence efficiency, particularly near the recovery threshold and under structured SORS measurements. Future work will further explore the interaction between structured modewise measurements and block-weighted Riemannian optimization, with the aim of developing more efficient and scalable methods for large-scale tensor recovery.

Appendix: Complete Convergence Proofs This appendix supplies the geometric and analytic estimates used in Section 4 and gives the complete proof of Theorem 4.3. (i) (i) For each mode i ∈ [d], we introduce the modewise orthogonal projectors P (i) and P U (i) , which Vl

(i)

project the mode-i fibers of any ambient tensor onto the column spaces of the current factor Vl and the (i) (i) true factor U (i) , respectively. They are equivalent to multiplying a tensor in mode i by Vl (Vl )T and (i) (i) T U (U ) . Lemma .1 (Perturbation of the mode subspaces). For every i ∈ [d], P

(i) (i)

(i)

− P U (i) ≤

∥Xl − T ∥F . σri (T(i) )

Proof. By the definition of the modewise projectors and the invariance of the Frobenius norm under (i) (i) (i) (i) (i) matricization, ∥P (i) −P U (i) ∥ = ∥Vl (Vl )T −U (i) (U (i) )T ∥2 . Since Vl and U (i) have orthonormal Vl

columns and span subspaces of the same dimension ri , the standard identity for equal-rank orthogonal (i) (i) projectors gives ∥Vl (Vl )T − U (i) (U (i) )T ∥2 = ∥P V (i)⊥ U (i) ∥2 . l

Let Ū (i) denote the Kronecker product of the true factor matrices excluding the i-th mode, defined (i) analogously to V̄l . Then T(i) = U (i) B(i) (Ū (i) )T . Since mulrank(T ) = r, B(i) has full row rank. Moreover, because U (i) and Ū (i) have orthonormal columns, T(i) and B(i) have the same nonzero singular † † values. Hence ∥B(i) ∥2 = 1/σri (T(i) ), and U (i) = T(i) Ū (i) B(i) . (i)

Since the column space of (Xl )(i) is contained in span(Vl ), P V (i)⊥ (Xl )(i) = 0. Therefore, l

† ∥P V (i)⊥ U (i) ∥2 = P V (i)⊥ (T − Xl )(i) Ū (i) B(i) l

† ≤ ∥P V (i)⊥ ∥2 ∥(T − Xl )(i) ∥F ∥Ū (i) ∥2 ∥B(i) ∥2 ≤ l

∥Xl − T ∥F . σri (T(i) )

For later use, we define the orthogonal projector associated with the row space of the i-th factor tangent block by P

(i)

(j)

Bl ,{Vl

}j̸=i

= Y(i) V̄l (Bl )†(i) (Bl )(i) (V̄l )T . (i)

Y (i)

(i)

(33)

Because (Bl )†(i) (Bl )(i) is the orthogonal projector onto the row space of (Bl )(i) and V̄l has orthonormal columns, (33) indeed defines an orthogonal projector. The closed-form expressions imply the following decomposition of the canonical tangent-space projector: d d Y X (i) (i) (i) P Sl = P (i) + P P (i)⊥ . (34) (j) (i)

i=1

Bl ,{Vl

}j̸=i

The projectors appearing in each product in (34) commute because the mode-i projector acts on the left (i) side of the mode-i matricization, whereas P acts on the corresponding right side. (j) Bl ,{Vl

}j̸=i

Lemma .2 (Second-order tangent-space mismatch). The orthogonal tangent-space projector satisfies 2d − 1 ∥Xl − T ∥2F . min1≤i≤d σri (T(i) )

∥(I − P Sl )T ∥F ≤ Proof. Expanding I =

(35)

(i) (i) i=1 (P V (i) + P V (i)⊥ ) and subtracting (34) decomposes I − P Sl into d

singleton terms (|Λ| = 1) and 2d − 1 − d higher-order terms (|Λ| ≥ 2):      d X Y (j) X Y Y (j) (j) (i)  P (j) − P (i) (j)  P (i)(i)⊥ +  I −P Sl = P (j)⊥   P (j)  P (i)⊥ . i=1

j̸=i

Bl ,{Vl

}j̸=i

Λ⊆[d] j∈Λ\{i} |Λ|≥2 |

}

:=E i

j ∈Λ /

:=Ee Λ

} (36)

where, for each higher-order term, i denotes an arbitrary index chosen from Λ. (i) (i) (i) (i) Using the identity P U (i) T = T , we have P (i)⊥ T = (P U (i) − P (i) )T for every i ∈ [d]. Vl

For the singleton terms, both projectors defining E i act as the identity on Xl , and hence E i Xl = 0. Q (i) (j) Moreover, the range of P is contained in the range of j̸=i P (j) . Therefore, E i is the (j) Bl ,{Vl

}j̸=i

difference of two nested orthogonal projectors and satisfies ∥E i ∥ ≤ 1. Since E i commutes with the mode-i projectors, E iP

(i) (i)⊥ Vl

(i)

T = (P U (i) − P

(i) (i) Vl

(i)

)E i T = (P U (i) − P

(i) (i)

)E i (T − Xl ).

Hence, Lemma .1 gives E iP

(i)

(i)⊥ T

≤ F

∥T − Xl ∥2F . σri (T(i) )

For each higher-order term, Ee Λ contains at least one mode-j orthogonal complement projector with j ̸= i, so that Ee Λ Xl = 0. Furthermore, since it is a product of commuting orthogonal projectors, ∥Ee Λ ∥ ≤ 1, and it commutes with the mode-i projectors. The same argument therefore yields (i) Ee Λ P (i)⊥ T

≤

∥T − Xl ∥2F . σri (T(i) )

There are d singleton terms and 2d − 1 − d higher-order terms in (36). Summing their Frobenius norms and uniformly bounding the denominators by min1≤i≤d σri (T(i) ) yields (35). Lemma .3 (Restricted tangent-space operator bounds). Assume that L satisfies TRIP(δ2r , 2r). Then ∥P Sl − P Sl L ∗ L P Sl ∥ ≤ δ2r ,

∥P Sl L ∗ L P Sl ∥ ≤ 1 + δ2r .

(37)

Proof. The first operator in (37) is self-adjoint and vanishes on Sl⊥ , so its operator norm is determined by its restriction to Sl . For any tensor Z ∈ Sl with ∥Z∥F = 1, we have P Sl Z = Z and mulrank(Z) ⪯ 2r, and hence the TRIP implies |∥Z∥2F − ∥L (Z)∥22 | ≤ δ2r . By the linearity of the inner product, the definition of the adjoint operator L ∗ , and the property P Sl Z = Z, we can expand the quadratic form as: ⟨(P Sl − P Sl L ∗ L P Sl ) Z, Z⟩F = ⟨P Sl Z, Z⟩F − ⟨P Sl L ∗ L P Sl Z, Z⟩F = ⟨Z, Z⟩F − ⟨L ∗ L Z, Z⟩F = ∥Z∥2F − ∥L (Z)∥22 ≤ δ2r . Thus, the first bound follows by taking the supremum of the absolute value over unit tensors in Sl . The second operator is positive semidefinite and also vanishes on Sl⊥ . By a similar Rayleigh quotient argument, the upper TRIP inequality guarantees ⟨P Sl L ∗ L P Sl Z, Z⟩F = ∥L (Z)∥22 ≤ 1 + δ2r for any unit tensor Z ∈ Sl , yielding ∥P Sl L ∗ L P Sl ∥ ≤ 1 + δ2r . 25

Lemma .4 (Rank and TRIP bound for the normal residual). The tensor (I − P Sl )T has multilinear rank at most 2r. Consequently, under TRIP(δ2r , 2r), it holds that ∥P Sl L ∗ L (I − P Sl )T ∥F ≤ (1 + δ2r )∥(I − P Sl )T ∥F . Proof. By the tangent projection formula (24), the mode-k column space of P Sl T is contained in the (k) sum of span(Vl ) and col(P V (k)⊥ T(k) ). Since col(T(k) ) = span(U (k) ), we have col(P V (k)⊥ T(k) ) ⊆ l

(k) (k) span(U (k) , Vl ). Hence col((P Sl T )(k) ) ⊆ span(U (k) , Vl ). Since the same inclusion holds for T(k) , (k) it follows that col(((I − P Sl )T )(k) ) ⊆ span(U (k) , Vl ). Therefore, mulrank((I − P Sl )T ) ⪯ 2r. For the cross-term bound, the self-adjointness of P Sl and the duality of the Frobenius norm give

∥P Sl L ∗ L (I − P Sl )T ∥F = sup |⟨L (I − P Sl )T , L (Z)⟩| ≤ (1 + δ2r ) ∥(I − P Sl )T ∥F , Z∈Sl ∥Z∥F =1

where the inequality follows from the Cauchy–Schwarz inequality and the upper TRIP bound applied separately to (I − P Sl )T and Z, both of which have multilinear rank at most 2r. We next give the complete proofs of the two adaptive perturbation estimates used in Section 4.2. Lemma .5 (Operator bound for the weight perturbation). Let ρl := max0≤k≤d |ωl,k − 1|. Then, for every ambient tensor Z, l (38) ∥(P ω l − P Sl )Z∥F ≤ ρl ∥P Sl Z∥F . In particular, l ∥(P ω l − P Sl )∥Sl →Sl ≤ ρl .

(0)

(d)

Proof. Since the ranges of Πl , . . . , Πl 2 l (P ω l − P Sl )Z F =

d X

(39)

are mutually orthogonal, for every ambient tensor Z we have

(k) (ωl,k − 1) ∥Πl Z∥2F ≤ ρ2l 2

k=0

d X

(k)

∥Πl Z∥2F = ρ2l ∥P Sl Z∥2F .

(40)

k=0

Taking square roots in (40) gives (38). For Z ∈ Sl , we have P Sl Z = Z, and hence taking the supremum over nonzero tensors in Sl gives (39). Lemma .6 (Adaptive-weight cross term). Assume that L satisfies TRIP(δ2r , 2r). Then ∗ l ∥(P ω l − P Sl )L L (Xl − T )∥F ≤ ρl (1 + δ2r ) ∥Xl − T ∥F +

2d − 1 ∥Xl − T ∥2F min1≤i≤d σri (T(i) )

(41)

Proof. Since Xl ∈ Sl , we have P Sl Xl = Xl and therefore Xl − T = P Sl (Xl − T ) − (I − P Sl )T . Applying Lemma .5 to L ∗ L (Xl − T ) and then using this decomposition gives ∗ l (P ω l − P Sl )L L (Xl − T ) F

≤ ρl ∥P Sl L ∗ L P Sl (Xl − T )∥F + ρl ∥P Sl L ∗ L (I − P Sl )T ∥F ≤ ρl (1 + δ2r ) (∥Xl − T ∥F + ∥(I − P Sl )T ∥F ) .

(42)

In the last inequality of (42), we used Lemma .3, Lemma .4, and ∥P Sl (Xl − T )∥F ≤ ∥Xl − T ∥F . Applying the second-order estimate (35) to the remaining normal component proves (41). Lemma .7 (Step-size bound for exact line search). Suppose ρl := max0≤k≤d |ωl,k − 1| < 1 and L satisfies TRIP(δ2r , 2r). Then the step size αl in (21) satisfies 1 1 ≤ αl ≤ . (1 + δ2r )(1 + ρl ) (1 − δ2r )(1 − ρl )

(43)

Consequently, |αl − 1| ≤

1 − 1. (1 − δ2r )(1 − ρl ) 26

(44)

(k)

Proof. Let al,k := ∥Πl Gl ∥F . The mutual orthogonality of the tangent blocks yields ⟨Gl , Zl ⟩F = Pd Pd 2 2 2 2 k=0 ωl,k al,k and ∥Zl ∥F = k=0 ωl,k al,k . Since 1 − ρl ≤ ωl,k ≤ 1 + ρl , it follows that (1 − ρl )⟨Gl , Zl ⟩F ≤ ∥Zl ∥2F ≤ (1 + ρl )⟨Gl , Zl ⟩F .

(45)

Moreover, Zl ∈ Sl implies mulrank(Zl ) ⪯ 2r. Hence TRIP(δ2r , 2r) gives (1 − δ2r )∥Zl ∥2F ≤ ∥L (Zl )∥22 ≤ (1 + δ2r )∥Zl ∥2F .

(46)

⟨Gl ,Zl ⟩F Combining (45) and (46) with αl = ∥L yields (43). (Zl )∥22 The upper endpoint in (43) gives the larger deviation from one, which proves (44). (k)

Lemma .8 (Uniform p control of the adaptive-weight deviation). Suppose that max0≤k≤d ∥Πl Gl ∥F ≤ Γ, and define SΓ := Γ2 + η 2 . Then ρl := max |ωl,k − 1| ≤ ρη (Γ) := 0≤k≤d

d(SΓ − η) . SΓ + dη

(47)

Proof. By the assumed block-gradient bound and the definition of sl,k , we have η ≤ sl,k ≤ SΓ for all k = 0, . . . , d. Hence, for each k, (d + 1)η (d + 1)SΓ ≤ ωl,k ≤ . η + dSΓ SΓ + dη d(SΓ −η) Γ −η) Therefore, ωl,k − 1 ≤ d(S SΓ +dη , 1 − ωl,k ≤ η+dSΓ . Since SΓ ≥ η and d ≥ 1, SΓ + dη ≤ η + dSΓ . Thus, Γ −η) |ωl,k − 1| ≤ d(S SΓ +dη = ρη (Γ). Taking the maximum over k = 0, . . . , d proves (47).

Under the TRIP(δ2r , 2r) assumption, whenever ∥Xl − T ∥F ≤ R, the block-gradient estimate gives (k)

∥Πl Gl ∥F =

sup (k)

|⟨Z, L ∗ L (Xl − T )⟩F | =

Z∈range(Πl ) ∥Z∥F =1

≤

|⟨L (Z), L (Xl − T )⟩|

sup (k)

Z∈range(Πl ) ∥Z∥F =1

∥L (Z)∥2 ∥L (Xl − T )∥2 ≤ (1 + δ2r )∥Xl − T ∥F .

sup

(48)

(k)

Z∈range(Πl ) ∥Z∥F =1

Applying Lemma .8 with Γ = (1 + δ2r )R yields ∥Xl − T ∥F ≤ R =⇒ ρl ≤ ρη ((1 + δ2r )R). If η is chosen such that ρη ((1 + δ2r )R) ≤ ρ⋆ , then ∥Xl − T ∥F ≤ R

=⇒

ρl ≤ ρ⋆ .

(49)

Complete proof of Theorem 4.3

Proof. Fix l ≥ 0 such that ∥Xl − T ∥F ≤ R. By (49), ρl ≤ ρ⋆ < 1. If Zl = 0, then P Sl Gl = 0. Lemmas .2 and .3 then force ∥Xl − T ∥F = 0 (since γR < 1), meaning Xl = T , and the theorem holds. Assume Zl ̸= 0. Since Zl ∈ Sl , mulrank(Zl ) ⪯ 2r. The lower TRIP bound guarantees ∥L (Zl )∥22 ≥ (1 − δ2r )∥Zl ∥2F > 0, ensuring αl is well-defined. Since Xl+1 = H r (Yl ), (6) and T ∈ Mr imply √ √ ∥Xl+1 − Yl ∥F ≤ d ∥Yl − P Mr (Yl )∥F ≤ d ∥Yl − T ∥F . Therefore, by the triangle inequality, √ ∥Xl+1 − T ∥F ≤ ∥Xl+1 − Yl ∥F + ∥Yl − T ∥F ≤ ( d + 1)∥Yl − T ∥F .

∗ l Substituting Yl = Xl − αl P ω l L L (Xl − T ) and utilizing Xl ∈ Sl , the error decomposes as √ ∗ l ∥Xl+1 − T ∥F ≤ ( d + 1) (Xl − T ) − αl P ω l L L (Xl − T ) F √ ∗ l = ( d + 1) (Xl − T ) − αl P Sl L ∗ L (Xl − T ) − αl (P ω l − P Sl )L L (Xl − T ) F √ ≤ ( d + 1) ∥(P Sl − αl P Sl L ∗ L P Sl ) (Xl − T )∥F + ∥(I − P Sl )(Xl − T )∥F | {z } | {z } I1

∗

+ αl ∥P Sl L L (I − P Sl )(Xl − T )∥F + αl | {z } |

∗ l (P ω l − P Sl )L L (Xl − T ) F

√ = ( d + 1)(I1 + I2 + I3 + I4 ).

{z I4

We next bound these four terms separately. • Bound of I1 : Writing P Sl −αl P Sl L ∗ L P Sl = (P Sl −P Sl L ∗ L P Sl )+(1−αl )P Sl L ∗ L P Sl , Lemma .3 and (44) yield δ2r + ρl − δ2r ρl 2δ2r + ρl − δ2r ρl I1 ≤ δ2r + (1 + δ2r ) ∥Xl − T ∥F = ∥Xl − T ∥F . (1 − δ2r )(1 − ρl ) (1 − δ2r )(1 − ρl ) • Bound of I2 : Applying Lemma .2 directly gives I2 = ∥(I − P Sl )T ∥F ≤

2d − 1 ∥Xl − T ∥2F . mini σri (T(i) )

• Bound of I3 : Lemmas .4 and .7 imply I3 ≤

2d − 1 1 + δ2r ∥Xl − T ∥2F . (1 − δ2r )(1 − ρl ) mini σri (T(i) )

• Bound of I4 : Combining Lemmas .6 and .7 leads to ρl (1 + δ2r ) 2d − 1 2 I4 ≤ ∥Xl − T ∥F + ∥Xl − T ∥F . (1 − δ2r )(1 − ρl ) mini σri (T(i) ) Combining bounds of I1 , I2 , I3 , and I4 yields " # √ (2d − 1)(1 + δ2r ρl ) 2( d + 1) ∥Xl+1 − T ∥F ≤ δ2r + ρl + ∥Xl − T ∥F ∥Xl − T ∥F . (50) (1 − δ2r )(1 − ρl ) mini σri (T(i) ) We proceed by induction. For l = 0, the assumption ∥X0 − T ∥F ≤ R and (49) ensure ρ0 ≤ ρ⋆ . Assume ∥Xl − T ∥F ≤ R, which implies ρl ≤ ρ⋆ . The coefficient multiplying ∥Xl − T ∥F in (50) is monotonically increasing with respect to both ∥Xl − T ∥F and ρl . Bounding them by R and ρ⋆ respectively gives: ∥Xl+1 − T ∥F ≤ γR ∥Xl − T ∥F . Since γR < 1, we obtain ∥Xl+1 − T ∥F ≤ R, completing the induction. Iterating this strict contraction l ∥X − T ∥ , and ρ ≤ ρ holds uniformly for all l ≥ 0. yields ∥Xl − T ∥F ≤ γR 0 ⋆ F l

}

Standard RGD initialization

The following result verifies the initialization condition used in Corollary 4.4 for the initialization prescribed in Algorithm 2. Lemma .9 (Initialization error). Let y = L (T ) and X0 = H r (L ∗ y). If L satisfies TRIP(δ2r , 2r), then (27) holds. Proof. For each i ∈ [d], let K (i) ∈ Rni ×ti , with ti ≤ 2ri , have orthonormal columns spanning the sum of span(U (i) ) and the subspace generated by the leading ri left singular vectors of (L ∗ y)(i) . Define the Q (i) orthogonal projector P K := di=1 P K (i) . By construction, P K T = T , and mulrank(P K Z) ⪯ 2r for any ambient tensor Z. According to the ordering property of the HOSVD [40], P K L ∗ y and L ∗ y share identical leading singular vectors and values, yielding X0 = H r (L ∗ y) = H r (P K L ∗ y). The quasi-optimality of the truncated HOSVD implies √ ∥X0 − P K L ∗ y∥F ≤ d∥T − P K L ∗ y∥F . Applying the triangle inequality, y = L (T ), and P K T = T , we obtain √ ∥X0 − T ∥F ≤ ∥X0 − P K L ∗ y∥F + ∥P K L ∗ y − T ∥F ≤ ( d + 1)∥T − P K L ∗ L T ∥F √ = ( d + 1)∥(P K − P K L ∗ L P K )T ∥F . Since P K − P K L ∗ L P K is self-adjoint and mulrank(P K Z) ⪯ 2r, the TRIP implies ∥P K − P K L ∗ L P K ∥ = sup |⟨(P K − P K L ∗ L P K )Z, Z⟩F | ∥Z∥F =1

= sup ∥Z∥F =1

Consequently,

∥P K Z∥2F − ∥L (P K Z)∥22 ≤ δ2r .

√ ∥X0 − T ∥F ≤ ( d + 1)δ2r ∥T ∥F .

References [1] Nan Li and Baoxin Li. Tensor completion for on-board compression of hyperspectral images. In 2010 IEEE International Conference on Image Processing, pages 517–520. IEEE, 2010. [2] Ji Liu, Przemyslaw Musialski, Peter Wonka, and Jieping Ye. Tensor completion for estimating missing values in visual data. IEEE Transactions on Pattern Analysis and Machine Intelligence, 35(1):208–220, 2013. [3] Zemin Zhang and Shuchin Aeron. Exact tensor completion using t-SVD. IEEE Transactions on Signal Processing, 65(6):1511–1526, 2017. [4] Kyle Gilman, Davoud Ataee Tarzanagh, and Laura Balzano. Grassmannian optimization for online tensor completion and tracking with the t-SVD. IEEE Transactions on Signal Processing, 70:2152– 2167, 2022. [5] Jie Lin, Li-Yuan Li, Xi-Le Zhao, and Han Yang. Robust multi-subspace tensor dictionary learning for multi-temporal image reconstruction. Inverse Problems, 42(7):075017, 2026. [6] Animashree Anandkumar, Rong Ge, Daniel Hsu, Sham M Kakade, and Matus Telgarsky. Tensor decompositions for learning latent variable models. Journal of Machine Learning Research, 15:2773–2832, 2014. 29

[7] Andrzej Cichocki, Danilo Mandic, Lieven De Lathauwer, Guoxu Zhou, Qibin Zhao, Cesar Caiafa, and Huy Anh Phan. Tensor decompositions for signal processing applications: From two-way to multiway component analysis. IEEE Signal Processing Magazine, 32(2):145–163, 2015. [8] Shipeng Zhang, Lizhi Wang, Lei Zhang, and Hua Huang. Learning tensor low-rank prior for hyperspectral image reconstruction. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 12006–12015, 2021. [9] David Gross, Yi-Kai Liu, Steven T Flammia, Stephen Becker, and Jens Eisert. Quantum state tomography via compressed sensing. Physical Review Letters, 105(15):150401, 2010. [10] Richard Kueng, Holger Rauhut, and Ulrich Terstiege. Low-rank matrix recovery from rank-one measurements. Applied and Computational Harmonic Analysis, 42(1):88–116, 2017. [11] Ledyard R Tucker. Some mathematical notes on three-mode factor analysis. Psychometrika, 31(3):279–311, 1966. [12] Simon Foucart. Sparse recovery algorithms: sufficient conditions in terms of restricted isometry constants. In Approximation Theory XIII: San Antonio 2010, pages 65–77. Springer, 2011. [13] Qun Mo and Song Li. New bounds on the restricted isometry constant δ2k . Applied and Computational Harmonic Analysis, 31(3):460–468, 2011. [14] Deanna Needell and Joel A Tropp. CoSaMP: Iterative signal recovery from incomplete and inaccurate samples. Applied and Computational Harmonic Analysis, 26(3):301–321, 2009. [15] Thomas Blumensath and Mike E. Davies. Iterative hard thresholding for compressed sensing. Applied and Computational Harmonic Analysis, 27(3):265–274, 2009. [16] Mark A Iwen, Deanna Needell, Elizaveta Rebrova, and Ali Zare. Lower memory oblivious (tensor) subspace embeddings with fewer random bits: modewise methods for least squares. SIAM Journal on Matrix Analysis and Applications, 42(1):376–416, 2021. [17] Ruhui Jin, Tamara G Kolda, and Rachel Ward. Faster Johnson–Lindenstrauss transforms via Kronecker products. Information and Inference: A Journal of the IMA, 10(4):1533–1562, 2021. [18] Cullen A Haselby, Mark A Iwen, Deanna Needell, Michael Perlmutter, and Elizaveta Rebrova. Modewise operators, the tensor restricted isometry property, and low-rank tensor recovery. Applied and Computational Harmonic Analysis, 66:161–192, 2023. [19] Stefan Bamberger, Felix Krahmer, and Rachel Ward. Johnson–Lindenstrauss embeddings with Kronecker structure. SIAM Journal on Matrix Analysis and Applications, 43(4):1806–1850, 2022. [20] Osman Asif Malik and Stephen Becker. Low-rank Tucker decomposition of large tensors using TensorSketch. In Advances in Neural Information Processing Systems, volume 31, pages 10096– 10106, 2018. [21] Holger Rauhut, Reinhold Schneider, and Željka Stojanac. Low-rank tensor recovery via iterative hard thresholding. Linear Algebra and its Applications, 523:220–262, 2017. [22] Silvia Gandy, Benjamin Recht, and Isao Yamada. Tensor completion and low-n-rank tensor recovery via convex optimization. Inverse Problems, 27(2):025010, 2011. [23] Bo Huang, Cun Mu, Donald Goldfarb, and John Wright. Provable low-rank tensor recovery. Optimization-Online, 4252(2):455–500, 2014.

[24] Cun Mu, Bo Huang, John Wright, and Donald Goldfarb. Square deal: Lower bounds and improved relaxations for tensor recovery. In Proceedings of the 31st International Conference on Machine Learning, volume 32 of Proceedings of Machine Learning Research, pages 73–81. PMLR, 2014. [25] Ming Yuan and Cun-Hui Zhang. On tensor completion via nuclear norm minimization. Foundations of Computational Mathematics, 16(4):1031–1068, 2016. [26] Ben-Zheng Li, Xi-Le Zhao, Hao Zhang, and Delin Chu. Importance-aware nonlocal tensor nuclear norm for high-dimensional image recovery. Inverse Problems, 42(1):015004, 2026. [27] Jian Lu, Lin Huang, Xiaoxia Liu, Ning Xie, Qingtang Jiang, and Yuru Zou. 3d poissonian image deblurring via patch-based tensor logarithmic schatten-p minimization. Inverse Problems, 40(6):065010, 2024. [28] Prateek Jain and Sewoong Oh. Provable tensor factorization with missing data. In Advances in Neural Information Processing Systems, volume 27, pages 1431–1439, 2014. [29] Dong Xia and Ming Yuan. On polynomial time methods for exact low-rank tensor completion. Foundations of Computational Mathematics, 19(6):1265–1313, 2019. [30] DONG XIA, MING YUAN, and CUN-HUI ZHANG. Statistically optimal and computationally efficient low-rank tensor completion from noisy entries. The Annals of Statistics, 49(1):76–99, 2021. [31] Daniel Kressner, Michael Steinlechner, and Bart Vandereycken. Low-rank tensor completion by riemannian optimization. BIT Numerical Mathematics, 54(2):447–468, 2014. [32] Jian-Feng Cai, Lizhang Miao, Yang Wang, and Yin Xian. Provable near-optimal low-multilinearrank tensor recovery. arXiv preprint arXiv:2007.08904, 2020. [33] Haifeng Wang, Jinchi Chen, and Ke Wei. Implicit regularization and entrywise convergence of Riemannian optimization for low Tucker-rank tensor completion. Journal of Machine Learning Research, 24(347):1–84, 2023. [34] Yuanwei Zhang, Fengmiao Bian, Xiaoqun Zhang, and Jian-Feng Cai. Preconditioned Riemannian gradient descent algorithm for low-multilinear-rank tensor completion. In Proceedings of the 42nd International Conference on Machine Learning, volume 267 of Proceedings of Machine Learning Research, pages 74490–74514. PMLR, 2025. [35] Fengmiao Bian, Jian-Feng Cai, and Rui Zhang. A preconditioned Riemannian gradient descent algorithm for low-rank matrix recovery. SIAM Journal on Matrix Analysis and Applications, 45(4):2075–2103, 2024. [36] Mohammad Hamed and Reshad Hosseini. Riemannian preconditioned coordinate descent for low multilinear rank approximation. SIAM Journal on Matrix Analysis and Applications, 45(2):1054– 1075, 2024. [37] Cullen Haselby, Mark A Iwen, Deanna Needell, Elizaveta Rebrova, and William Swartworth. Fast and low-memory compressive sensing algorithms for low Tucker-rank tensor approximation from streamed measurements. Numerical Algorithms, 101(4):2567–2629, 2026. [38] Bamdev Mishra and Rodolphe Sepulchre. Riemannian preconditioning. SIAM Journal on Optimization, 26(1):635–660, 2016. [39] Nick Vannieuwenhoven, Raf Vandebril, and Karl Meerbergen. A new truncation strategy for the higher-order singular value decomposition. SIAM Journal on Scientific Computing, 34(2):A1027– A1052, 2012. 31

[40] Lieven De Lathauwer, Bart De Moor, and Joos Vandewalle. A multilinear singular value decomposition. SIAM Journal on Matrix Analysis and Applications, 21(4):1253–1278, 2000. [41] Christopher J Hillar and Lek-Heng Lim. Most tensor problems are NP-hard. Journal of the ACM (JACM), 60(6):1–39, 2013.

Record · ID 1028703 · SHA-256 9a12616384a810ac
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.