{"id":"136bc82c-6e82-4af6-b631-24f1a988d06c","arxiv_id":"2506.01658","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":4,"one_line_summary":"A two-sided CP-decomposed structure for tensor autoregression regresses response features on covariate features, matching Tucker-level interpretability at CP-level parameter counts, plus a low-rank-plus-sparse variant with non-asymptotic error bounds.","lead":"The paper proposes an autoregressive model for tensor-valued time series whose coefficient tensor is built from two CP-decomposed loading sets, one compressing the response into features and the other compressing past values into covariate features, which are then regressed on each other.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Non-asymptotic bounds rely on covering-number estimates that require a uniform bound on representation norms of G, which Assumptions 1-7 do not supply.","rationale":"The reader's weakest assumption identifies the same point I reach: the uniform conditioning/representation-norm bound is absent from Assumptions 1–7. I agree with the CONDITIONAL verdict because the contribution is real—the two-sided CP parameterization, the linear-in-order parameter count, and the ALS algorithm are valuable—and the proof gap is likely repairable by adding a lower bound on the smallest singular values of Λ_y and Λ_x (or an equivalent upper bound on representation norms of G). This would not change the rates but would make the theorems valid. The simulations and ENSO analysis are supportive but not decisive for the theoretical claim; the absence of error bars on forecast comparisons and the flat FPR in Table 1 are secondary concerns. The most load-bearing issue remains the covering-number lemma, and since the reader already flagged it, my recommendation is UNCHANGED.","tokens_in":42372,"tokens_out":13908,"duration_ms":160338,"concrete_test":"Verify Lemma 1(b) in the minimal tensor case n=3, q=(2,2,2), P=1, Ry=2, Rx=1. Construct unit CP vectors a_ε, b_ε ∈ R^8 with |⟨a_ε,b_ε⟩|=1−ε (e.g., by taking the three factor pairs nearly parallel), set Δ_ε = (a_ε − b_ε)/‖a_ε − b_ε‖. If Lemma 1(b) is correct, Δ_ε ∈ Ξ' must be within ε of the net (B.1), whose G-ball has radius √2 g·√min(2,1) = √2 g. Compute (symbolically or via a rank-2 CP decomposition solver) the minimal ‖G‖_F over all representations Δ_ε = Λ_y G Λ_x^T with Λ_y ∈ S((2,2,2),2), Λ_x ∈ S((2,2,2),1). If this minimum diverges as ε→0, the covering bound in Lemma 1(b) fails as stated without an extra conditioning assumption.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central claim—that the OLS estimator in (2.10) satisfies the non-asymptotic rate of Theorem 2, and the low-rank-plus-sparse estimator satisfies Theorem 4—is proved by generic chaining over the difference set Ξ' = {Δ = A1 − A2 : Ai ∈ Θ, ‖Δ‖_F ≤ 1}. The chaining argument needs the covering-number bounds of Lemma 1(b)–(d), whose proofs cover the middle factor G inside a ball of radius proportional to g·√min(Ry, P Rx). That radius is justified by Assumption 3 (‖G‖_op ≤ g). But Assumption 3 constrains only the true coefficient's G; the parameter space Θ in (2.10) contains every A of the form Λ_y G (I_P ⊗ Λ_x^T) with Λ_y,Λ_x ∈ S and G arbitrary. Because the columns of Λ_y and Λ_x are only unit-norm, not uniformly well-conditioned, elements of Θ can have representations with ‖G‖_F arbitrarily large even when ‖A‖_F ≤ 1: as two Khatri–Rao columns become nearly collinear, one must use large G entries to represent a small difference. The covering net in (B.1) then misses these elements, and the Gaussian-width estimate (B.14) is not justified. The reader is right that a lower bound on σ_min(Λ_y), σ_min(Λ_x), or an explicit uniform bound on the minimal representation norm of G, is missing from Assumptions 1–7. This is a genuine gap in the proofs of Theorems 1–4, though it is likely patchable by adding such an assumption.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"This paper proposes a new low-rank structure for the coefficient tensor in tensor-on-tensor regression and tensor autoregression. The structure writes the sequential matricization of the coefficient as Λ_y G Λ_xᵀ (with a Kronecker factor for lags), so that response and covariate features are extracted by CP-type loadings and linked by a small coefficient matrix G. This yields a parameter count linear in tensor order and a supervised-factor interpretation. The paper also adds a sparse component (low-rank plus sparse) with an identifiability radius α_L, gives OLS and penalized estimators, an alternating least-squares algorithm, nonasymptotic bounds (Theorems 1–4), simulations, and an ENSO application.","tokens_in":42610,"tokens_out":16969,"duration_ms":171265,"significance":"If the rate results are correct, the paper's contribution is substantial: it provides an interpretable alternative to Tucker-based tensor AR with dramatically smaller parameter count, and the low-rank-plus-sparse extension addresses misspecification. The paper is careful to state the sample-size-independent unidentifiability error in Theorem 4, and the supplement contains detailed proofs. The main obstacle is the missing uniform conditioning assumption in the covering-number arguments; this is a fixable gap but is load-bearing for the stated theorems.","major_comments":[{"comment":"The proofs of the covering-number bounds bound the middle factor by 'By Assumption 3, ∥G_1∥_F ≤ ... and same for G_2' for arbitrary A_1, A_2 in the parameter space. Assumption 3, however, bounds only the G of the true coefficient tensor, not all elements of Θ. Since S(q,R) only constrains columns to be unit-norm and does not bound singular values of Λ_y or Λ_x away from zero, there exist elements of Θ with ∥A∥_F ≤ 1 for which every representation requires arbitrarily large ∥G∥_F (nearly collinear Khatri–Rao columns). The nets in (B.1) and the Gaussian-width estimate (B.14) therefore do not cover the difference set Ξ', and the RSC/deviation lemmas are not justified. Please add a uniform conditioning assumption (e.g., σ_min(Λ_y), σ_min(Λ_x) ≥ c > 0 for all admissible loadings, or an explicit bound on the minimal representation norm of G) and verify that Lemmas 1–7 and Theorems 1–4 follow.","section":"Section B, Lemma 1(b)–(d), used in Theorems 1–4"},{"comment":"The decomposition uses (Λ_yᵀΛ_y)⁻¹, so Λ_y must have full column rank for the model to be well defined. The 'non-degeneracy' condition stated after Eq. (2.3) is not incorporated into the formal parameter space S(q,R) or into Assumptions 1–7; moreover, after the reparameterization G := (Λ_yᵀΛ_y)⁻¹G, the quantity bounded in Assumption 3 is not the original G. Please make the full-rank/quantitative conditioning requirement a formal assumption; this also helps resolve the issue in the previous comment.","section":"Eq. (2.4) and Section 2.3"}],"minor_comments":[{"comment":"The proof defines g′ with an extra factor P relative to the statement of the lemma; please reconcile the two definitions so the constants in the covering bound are consistent.","section":"Section B, Lemma 1(b)"},{"comment":"Because c_x and C_x are allowed to depend on dimensions and to shrink or diverge, the phrases 'bounded away from zero and infinity' in the discussion after Theorems 1–3 need a precise convention; otherwise the displayed rates hide dimension-dependent constants.","section":"Assumption 1 and Theorems 1–3"},{"comment":"The statement lists only Assumptions 3–5 in its hypotheses; include Assumptions 6 and 7 explicitly for completeness.","section":"Theorem 4"},{"comment":"The update of A_L via Algorithm 1 does not enforce the constraint ∥A_L∥_∞ ≤ α_L/(P Q²) from (2.11); clarify whether a final trimming step is applied or whether the constrained subproblem is solved.","section":"Algorithm 2"},{"comment":"The text abbreviates the error as MAFE while the table header reads MSAE; make the abbreviation consistent.","section":"Section 5, Table 2"},{"comment":"The notation ∥δ∥_2 = ∥∆^j∥_2 refers to an undefined ∆^j; define this row quantity explicitly.","section":"Section B, Lemma 1(d)"}],"recommendation":"major_revision","confidential_remarks":"The manuscript fits the journal's scope, and the main gap identified in the report is patchable by adding a uniform conditioning assumption. I do not see grounds for rejection, and I have no concerns about citation behavior."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Dear colleague,\n\nThe one thing you should know: this paper proposes a genuinely new low-dimensional structure for tensor-on-tensor regression and tensor autoregression — a two-sided CP decomposition with a pseudo-inverse normalization — and it deserves a serious referee, but the non-asymptotic proofs as written have a real gap that needs fixing.\n\nWhat is actually new: the structure [A]_n = Λ_y (Λ_y^T Λ_y)^{-1} G Λ_x^T makes the response-feature regression exact, giving a supervised factor interpretation at CP-level parameter counts. It reduces to CP tensor-on-tensor regression when R_y = R_x and G is diagonal, and to a constrained reduced-rank VAR otherwise. Parameter count is linear in tensor order, which matters for order-6 coefficient tensors like the ENSO example. That is a real middle ground and I have not seen it in the cited prior work.\n\nWhat the paper does well: the ALS algorithm is clearly laid out, the simulation study checks the claimed convergence rates, and the ENSO application yields interpretable spatial loadings. The paper is honest about limitations: it discloses the T ≳ sum q_i requirement and the sample-size-independent unidentifiability error in the low-rank-plus-sparse bound.\n\nNow the soft spots, in proportion. The main one is real: the covering-number Lemma 1(b)–(d) bounds the covering number of the difference set Ξ' by covering Λ_y, Λ_x, and G in a ball of radius proportional to g√min(R_y, P R_x). But the parameter space Θ only requires unit-norm columns for the loadings and no lower bound on their conditioning. When Khatri–Rao columns are nearly collinear, the minimal-norm G for a small difference can become arbitrarily large, so the nets miss elements of Ξ'. Assumptions 1–7 do not rule this out. This is a genuine gap in the proofs as written, though it looks patchable by adding a uniform lower bound on σ_min(Λ_y) and σ_min(Λ_x), or an explicit bound on representation norms of G. The rates themselves are plausible and standard in form.\n\nMinor concerns: no error bars on the forecast comparison, the margin over Lasso is thin, and the low-rank-plus-sparse model inherits ranks from the low-rank model rather than selecting its own. These are secondary.\n\nBottom line: this paper is worth refereeing. The model class is useful, the interpretation is new, and the theoretical gap is identifiable and likely fixable. I would accept it for peer review with the expectation of a revision that patches the conditioning assumption and tightens the empirical presentation.","headline":"A genuinely new two-sided CP structure for tensor autoregression, worth refereeing despite a real but likely patchable gap in the covering-number proofs.","tokens_in":43290,"tokens_out":3242,"would_cite":true,"duration_ms":33580,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["62M10","62J07","15A69"],"pacs":[],"model":"deepseek-v4-flash","headline":"The paper claims that a CP-based factorization of the coefficient tensor makes tensor autoregression interpretable as supervised factor regression while keeping parameters linear in tensor order.","keywords":["tensor-valued time series","autoregression","CP decomposition","supervised factor model","low-rank plus sparse","nonasymptotic bounds","alternating least squares","El Niño-Southern Oscillation"],"falsifier":"Simulate a stationary tensor AR(1) model with moderate dimensions, choose loading matrices whose columns are nearly collinear (angle around one degree), and compare the OLS estimation error against the predicted $\\sqrt{d_{\\mathrm{AR}} \\log(\\cdot)/T}$ rate as the angle decreases; if the error grows substantially, the rate is carried by an unstated uniform representation-norm condition rather than the assumptions listed.","tokens_in":41999,"feed_emoji":"📊","tokens_out":6995,"duration_ms":70949,"temperature":0.7,"pith_summary":"Tensor-valued time series are hard to model because the coefficient tensor of even an order-one autoregression has twice the data dimension, making direct estimation hopeless. The paper proposes a middle path between Tucker and CP decompositions: it factors the coefficient tensor so that response features extracted by one set of CP loadings are regressed on covariate features extracted by another set. This supervised-factor structure keeps the parameter count linear in tensor order, specifically $P R_y R_x + (R_y + R_x) \\sum_i q_i$ for an AR($P$) model, and the paper derives nonasymptotic ordinary-least-squares error bounds that match this complexity. A low-rank-plus-sparse extension is added to absorb heterogeneous or spiky signals that an exact low-rank assumption would miss. The motivating application is forecasting ocean temperature anomalies related to El Niño, where the fitted loadings recover the El Niño basin.","feed_headline":"Tensor AR gets interpretable factors with linear parameter growth","feed_subtitle":"The coefficient tensor is factored so response features regress on covariate features with linear parameter growth.","key_machinery":"The load-bearing object is the decomposition in equations (2.4)-(2.5), $[A]_n = \\Lambda_y (\\Lambda_y^\\top \\Lambda_y)^{-1} G \\Lambda_x^\\top$, which rewrites tensor-on-tensor regression as a regression of extracted response features on extracted covariate features. The loading matrices are Khatri-Rao products of unit-norm factor columns, giving each column the form of a vectorized CP rank-one tensor, and in the autoregressive version the same $\\Lambda_y$ and $\\Lambda_x$ are shared across lags with only $G_k$ varying. The proofs use covering-number bounds for the difference set of two such coefficient tensors together with restricted strong convexity of the empirical loss, and the optimization is an alternating least squares routine in which each block update reduces to a closed-form linear regression.","core_discovery":"The central claim is that the coefficient tensor of a tensor-on-tensor regression or tensor autoregression can be parameterized, via sequential matricization, as $[A]_n = \\Lambda_y (\\Lambda_y^\\top \\Lambda_y)^{-1} G \\Lambda_x^\\top$, where $\\Lambda_y$ and $\\Lambda_x$ are Khatri-Rao products of unit-norm loading matrices. Multiplying the regression equation by $\\Lambda_y^\\top$ yields $f_t^{\\mathrm{response}} = G f_t^{\\mathrm{covariate}} + e_t^{\\mathrm{CP}}$, so the model literally regresses response features on covariate features. The paper argues that this CP-based low-rank structure combines the interpretability of supervised factor models with the parsimony of CP decomposition, and that the OLS estimator over this parameter space attains Frobenius-error bounds of order $\\sqrt{d_{\\mathrm{AR}} d_c / (T-P)}$ under stationarity and sub-Gaussian noise. For the low-rank-plus-sparse variant, the error splits into an estimation term and a sample-size-independent unidentifiability term, the latter governed by the radius $\\alpha_L$.","pith_inferences":["The same factor-regression parametrization applies to any tensor-on-tensor regression, so the interpretability-efficiency trade-off it resolves is not specific to time series; one could test it directly on cross-sectional imaging or marketing data.","Sharing $\\Lambda_y$ and $\\Lambda_x$ across lags presumes a single linear feature extractor for all temporal dynamics; regime-switching or time-varying dynamics would likely require either multiple factor sets or an extension where loadings evolve.","The theorem's unidentifiability term suggests that any practical use of the low-rank-plus-sparse model should report the chosen $\\alpha_L$ together with the fitted sparse tensor, since the same likelihood can be explained by shifting mass between the two components.","A natural stress test of the theory is to let the loading columns become progressively closer to collinear and check whether the empirical error degrades as the covering-number argument predicts; the paper's assumptions do not rule out such a degradation."],"forward_implications":["For an AR($P$) model the parameter count is $d_{\\mathrm{AR}} = P R_y R_x + (R_y + R_x) \\sum_i q_i$, so the number of parameters grows linearly with tensor order instead of exponentially as in Tucker-based tensor AR.","The OLS estimator satisfies $\\|\\hat{A} - A^*\\|_F \\lesssim (\\kappa^2 \\kappa'_{U,B} / \\kappa_{L,A}) \\sqrt{d_{\\mathrm{AR}} d_c / (T-P)}$ with high probability, making the rate explicit in the sample size, ranks, and dimensions.","The low-rank-plus-sparse estimator separates low-rank collaborative structure from sparse idiosyncratic structure, with squared error bounded by $(d_{\\mathrm{AR}} d_c + s \\log(P Q^2))/(T-P)$ plus an unidentifiability term $s \\alpha_L^2/(P Q^2)$ that does not vanish with sample size.","The fitted response loadings from the ENSO example concentrate weight on the El Niño basin, giving the model a concrete spatial interpretation that pure CP or Lasso estimates lack."],"supporting_citations":[{"why":"Establishes the Tucker-based supervised factor interpretation for tensor autoregression that this paper builds on and contrasts with.","marker":"(Wang et al., 2024)"},{"why":"Provides the CP-based scalar-on-tensor regression framework and the alternating least squares estimation strategy adapted here.","marker":"(Zhou et al., 2013)"},{"why":"Supplies the spectral-measure and restricted-strong-convexity machinery for nonasymptotic analysis of high-dimensional VAR models.","marker":"(Basu and Michailidis, 2015)"},{"why":"Reviews CP and Tucker decompositions and defines the Khatri-Rao product conventions used in the parametrization.","marker":"(Kolda and Bader, 2009)"},{"why":"Provides Gaussian-width, $\\gamma_\\alpha$-functional, and Hanson-Wright tools used in the covering-number and deviation bounds.","marker":"(Vershynin, 2018)"},{"why":"Introduces the low-rank-plus-sparse decomposition and the radius-of-nonidentifiability idea that the paper adapts for the tensor AR setting.","marker":"(Basu et al., 2019)"}],"fun_headline_variants":["Tensor AR gets interpretable CP factors with linear growth","CP-based tensor AR: response features regress on covariate features","Low-rank plus sparse tensor AR for robust interpretability","Efficient tensor AR with supervised factor interpretation"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The theorems rely on every coefficient tensor in the parameter space having a representation whose inner matrix $G$ is bounded in norm; the paper only assumes non-degeneracy of the loadings, which excludes exact linear dependence but permits nearly collinear columns where that norm can diverge.","fun_headline_variants_meta":{"raw":{"variants":["Tensor AR gets interpretable CP factors with linear growth","CP-based tensor AR: response features regress on covariate features","Low-rank plus sparse tensor AR for robust interpretability","Efficient tensor AR with supervised factor interpretation"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.00019,"raw_usage":{"total_tokens":1352,"prompt_tokens":968,"completion_tokens":384,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":584,"completion_tokens_details":{"reasoning_tokens":320}},"tokens_in":584,"tokens_out":384,"duration_ms":4472,"temperature":1.0,"reasoning_tokens":320,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-07T11:40:11.829835+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Simulate a stationary tensor AR(1) model with moderate dimensions, choose loading matrices whose columns are nearly collinear (angle around one degree), and compare the OLS estimation error against the predicted $\\sqrt{d_{\\mathrm{AR}} \\log(\\cdot)/T}$ rate as the angle decreases; if the error grows substantially, the rate is carried by an unstated uniform representation-norm condition rather than the assumptions listed.","supporting_citations":[{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Introduces the low-rank-plus-sparse decomposition and the radius-of-nonidentifiability idea that the paper adapts for the tensor AR setting."}],"review_version":1}