{"id":"9bb95853-6afc-4046-a327-3b394fea2458","arxiv_id":"2411.18210","paper_version":1,"verdict":"REJECT","confidence":"MODERATE","novelty_score":2.0,"correctness_risk":"high","formal_verification":"none","parameter_count":2,"one_line_summary":"The paper derives tensor versions of rational Arnoldi and Lanczos algorithms for model reduction and tensor Lyapunov equations, but they are equivalent to existing matrix methods after unfolding.","lead":"This paper applies rational Krylov subspace methods, reformulated with tensor products, to reduce the order of multidimensional linear dynamical systems and to solve tensor Lyapunov equations. The authors report numerical tests, but the core result reduces to known matrix methods via tensor unfolding and key proofs contain errors.","discovery_kind":"extension","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Theorem 6.1's residual derivation assumes V_m W_m^T B = B, which fails under Algorithm 4's rational Lanczos initialization, invalidating the Lyapunov residual bound.","rationale":"The reader's weakest_assumption pinpoints the exact invalid step in Theorem 6.1: the assumption V_m W_m^T B = B. My stress-test confirms this is load-bearing. In Algorithm 4, the initial block V1 is a QR factor of (A - sigma1 I)^{-1} B, not B, so range(B) is generally not in the rational Krylov subspace spanned by V_m. The cancellation in the proof of Theorem 6.1 is therefore unjustified, and the derived reduced Lyapunov equation (6.5) does not follow from the stated Galerkin condition. This is a genuine mathematical flaw, not merely a disagreement with a community convention. It attacks the core algorithmic claim for solving tensor Lyapunov equations, which is a stated application of the paper. The novelty and reproducibility concerns raised by the reader are secondary: even if the methods were new and code were provided, this incorrect residual derivation would leave the central application unsupported. I therefore see no reason to change the reader's REJECT verdict; the concern is severe and would require a corrected initialization (e.g., including B in V1) or a different proof strategy to be resolved.","tokens_in":20545,"tokens_out":5137,"duration_ms":40931,"concrete_test":"Construct a small random test problem, e.g., A a 2x2x2x2 tensor and B a 2x2x1x1 tensor. Run Algorithm 4 for m=2 with distinct shifts sigma1 and sigma2, and compute V2 and W2. Check whether ||V2 * W2^T * B - B||_F is machine zero. For generic A and B it will be nonzero, demonstrating that the identity used in Theorem 6.1 fails. Then solve the reduced Lyapunov equation (6.5) for Y2, form X2 = V2 * Y2 * V2^T, and compute the residual R2 = A * X2 + X2 * A^T + B * B^T. Verify whether the claimed bound ||R2|| <= 2||Gamma_2,A X2 V2^T|| actually holds. If it does not, Theorem 6.1 is invalidated.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central correctness flaw identified by the reader is in Section 6.1, specifically Theorem 6.1 and the derivation of (6.5). Algorithm 4 (TRBL) initializes with S0 = (A - sigma1 I)^{-1} * B and takes a QR decomposition S0 = V1 * H1,0 (steps 3-4). Therefore V1 spans range((A - sigma1 I)^{-1} B), not range(B). Since subsequent V_k are generated from the rational Krylov subspace built from this V1, generically range(B) is not contained in range(V_m). The proof of Theorem 6.1 then cancels B*B^T by invoking 'Vm * Wm^T * B = B'. This identity is true for the standard Lanczos process that starts with B itself, but it does not hold for the rational Lanczos process unless B is in the range of (A - sigma1 I)^{-1} B, which is not true for generic A and B. Consequently, after substituting (3.7), the residual Rm contains an uncancelled term B*B^T - Vm*Wm^T*B*B^T*Wm*Vm^T, which is generically nonzero. Thus the low-dimensional Lyapunov equation (6.5) is not the Galerkin projection of (6.2), and the bound ||Rm|| <= 2||Gamma_m,A Xm Vm^T|| in (6.6) is unsupported. This undermines the paper's claimed application of the tensor rational Lanczos method to large-scale Lyapunov tensor equations, a central advertised contribution.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"This manuscript extends rational block Krylov subspace methods to tensors via the Einstein product, targeting multidimensional linear time-invariant (MLTI) systems. It presents a tensor rational block Arnoldi algorithm (Algorithm 3) and a tensor rational block Lanczos algorithm (Algorithm 4), together with projection-based reduced models, transfer-function error estimates (Theorems 3.1 and 3.2), an adaptive shift-selection procedure (Section 4), and an application to continuous-time Lyapunov tensor equations with a residual bound (Theorem 6.1). Numerical experiments in Section 7 compare the proposed methods with classical tensor Krylov variants and report residual norms, iteration counts, and runtimes. The paper concludes that the adaptive tensor rational methods are effective for model reduction and for Lyapunov tensor equations.","tokens_in":20792,"tokens_out":18817,"duration_ms":153246,"significance":"The paper addresses a timely problem: extending rational Krylov techniques to tensor-structured dynamical systems, where the Einstein product formulation avoids explicit matricization. The tensor rational block Arnoldi construction in Section 3.1 and the associated transfer-function error bound in Theorem 3.1 are structurally plausible and could be useful for MLTI model reduction. The paper also proposes a concrete Lanczos-type algorithm and reports experiments, and it attempts to provide a cheap residual estimator for Lyapunov tensor equations. However, the analysis of the rational Lanczos variant contains a load-bearing projection error that invalidates several theoretical claims, and several implementation details in Section 7 are not reproducible as written. If the Lanczos error estimates and the Lyapunov residual bound can be repaired, the TRBA portion and the overall framework would be a worthwhile contribution; in its present form the paper cannot be recommended.","major_comments":[{"comment":"The proof of Theorem 6.1 cancels the term B * B^T by invoking the identity Vm * Wm^T * B = B. Under Algorithm 4, the initial block is S0 = (A - sigma1 I)^{-1} * B and V1 is obtained from the QR factorization of S0 (steps 3-4), so B is not generically in Range(V_m). The identity is therefore false for this algorithm, and the residual expansion in the proof is missing the term B * B^T - Vm * Bm * Bm^T * Vm^T. Consequently the bound ||Rm|| <= 2 ||Gamma_m,A * Xm * Vm^T|| in (6.6) is unsupported, and the residual-based stopping criterion in Algorithm 5 is invalid. I note that equation (6.5) does follow directly from the Galerkin condition (6.4) without this identity, so the reduced Lyapunov equation itself is not in question; only the residual estimator and its use in the numerical experiments are.","section":"Section 6.1, Theorem 6.1"},{"comment":"The residuals RB(s) and RC(s) in (3.11)-(3.12) are derived by cancelling B - Vm * Bm, which requires B = Vm * Wm^T * B. Because Algorithm 4 starts with (A - sigma1 I)^{-1} * B rather than with B itself, this cancellation is not valid; the residual contains the additional term B - Vm * Wm^T * B (and its transpose analogue), which is nonzero in general. As a result, the simplified error expression (3.13), Theorem 3.3, and the adaptive shift selections (4.1)-(4.2) built on these residuals are not established for the rational Lanczos reduced model. The authors should either change Algorithm 4 to start with B and C^T, as is standard in the matrix rational Lanczos literature (e.g., Ref. [4]), or re-derive the residual formulas without this cancellation.","section":"Section 3.2, equations (3.11)-(3.12) and Theorem 3.2"},{"comment":"The text says 'We use Algorithm 3 to solve the two continuous Lyapunov equations (5.1)' in both Example 2 and Example 3. Algorithm 3 is the tensor rational block Arnoldi algorithm for model reduction and does not solve Lyapunov equations; no variant of Algorithm 3 for Lyapunov equations is described in the paper. If Algorithm 5 or Algorithm 6 was intended, the reference must be corrected. In addition, the convergence criterion in those algorithms relies on Theorem 6.1, whose residual bound is invalid per the comment above, so the residual values in Table 7.1 cannot be interpreted as evidence of convergence.","section":"Section 7, Examples 2 and 3"},{"comment":"The adaptive shift selection in (4.1)-(4.2) maximizes over a set S that is never defined. The paper should specify whether S is a finite grid, a subset of the imaginary axis, or a frequency interval, and how it is chosen and updated during the iteration; without this specification the adaptive procedure is not reproducible.","section":"Section 4"}],"minor_comments":[{"comment":"The title contains a grammatical error ('A rational Krylov methods'); it should be 'Rational Krylov methods' or 'A rational Krylov method'. The keyword 'Model recuction' on the first page is a typo for 'Model reduction'.","section":"Title and Abstract"},{"comment":"The heading '4. pole selection' should be capitalized: 'Pole selection'.","section":"Section 4"},{"comment":"Step 5 of Algorithm 5 instructs the reader to solve the continuous-time Lyapunov equation (6.5) with the MATLAB function dlyap, which solves the discrete-time Lyapunov equation; the continuous-time solver lyap (or another appropriate routine) should be used.","section":"Algorithm 5"},{"comment":"The example numbering is inconsistent: 'Example 1.1' appears twice, while the second experiment is later referenced as 'Example 1.2' in the captions of Figures 7.3 and 7.4. Also, the captions of Figures 7.2 and 7.4 mix the notations F and H for the transfer function.","section":"Section 7"},{"comment":"The definition of the tensor rational Krylov subspace in (3.1) lists only (A - sigma1 I)^{-1} * B and the full product of shifted inverses applied to B, but Algorithm 3 generates a basis that also contains B itself (via the initial QR of B) and intermediate products; the definition should be aligned with the algorithm.","section":"Section 3.1"},{"comment":"Theorem 3.2 is stated without proof, and the text then moves to a different residual formulation; a proof or a precise reference to the matrix analogue should be provided.","section":"Theorem 3.2"},{"comment":"In (5.3), the dimensions 'r2' and 'N2' are not defined; presumably r and J1*J2 are intended. The phrase 'weakly symmetric positive-definite square tensor' should be defined or cited precisely.","section":"Section 5"},{"comment":"For reproducibility, the random tensors used in the experiments should be generated with a fixed seed, and the tensorization of the matrix A should be described precisely.","section":"Section 7"},{"comment":"In the proof of Theorem 6.1, 'As Ym is a symmetric matrix' should read 'symmetric tensor'.","section":"Theorem 6.1 proof"},{"comment":"There are many typographical errors ('F orj', 'end F or', 'in this cas', 'Lya punov matrix equations' in Section 2.2); a careful proofread is needed.","section":"Throughout"}],"recommendation":"major_revision","confidential_remarks":"The paper fits the scope of math.NA, but the novelty relative to the existing tensor Krylov literature (e.g., Refs. [10] and [19]) is incremental, and the rational Lanczos analysis is currently not reliable. If the authors resubmit, the editor should require that the Lanczos starting-block issue be resolved, that the numerical experiments refer to the correct algorithms and use reliable residual computations, and that the notation and presentation be cleaned up. The referee does not see grounds for rejecting the TRBA portion outright, but the paper as a whole is not acceptable in its present form."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"The reader's rejection is on the right track, and I want to add one thing: the flaw is real and it runs deeper than Theorem 6.1. The identity V_m * W_m^T * B = B is not just used in the Lyapunov residual proof; it also underlies the rational Lanczos error expressions in Section 3.2 (eq. 3.11). Since Algorithm 4 initializes S0 = (A - σ1 I)^{-1} * B and then takes a QR of S0 to get V1, the range of V_m is a rational Krylov space built from a shifted solve, not from B. Unless σ1 = ∞ (i.e., an ordinary block Lanczos start), B does not generally lie in that range. So the Galerkin condition (6.4) does not give the reduced Lyapunov equation (6.5), and the residual bound (6.6) is unsupported.\n\nCredit where due: Section 3.1 is honest that the tensor results follow from the unfolding isomorphism, and the adaptive shift selection is a reasonable heuristic, though it follows prior matrix work. The paper could be useful if the proof were fixed or the algorithm modified to ensure B lies in the initial subspace.\n\nThat said, the novelty is thin: through the isomorphism Ψ, the methods are the matrix ones on an unfolded system, so the practical payoff must be argued with experiments, and that is exactly where the paper is weak. The numerical experiments give no code or data, the search set S for the adaptive shifts is never defined, and there is no comparison against simply unfolding the tensors and applying standard matrix rational Krylov—the natural baseline given the isomorphism. Table 7.1 shows TRBL is about 10x slower than TCBL, so the claimed effectiveness is not obvious. Also, Algorithm 5 calls dlyap for a continuous Lyapunov equation; that is a discrete-time solver.\n\nMy take: this deserves a serious referee, but the current form should not be published. The false identity is load-bearing, and the experiments are not yet convincing. I would not cite it until the derivation is corrected and the experiments are reproducible. For a reading group, it might be a useful case study in how biorthogonality conditions can be misapplied to rational Krylov bases, but otherwise it is a niche contribution with a fixable flaw.","headline":"The tensor rational Lanczos sections rest on an identity that fails for the algorithm's own starting block; the paper is honest about its lineage but the Lyapunov result is currently unsupported.","tokens_in":21393,"tokens_out":6119,"would_cite":false,"duration_ms":52889,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["15A69","93B11"],"pacs":[],"model":"deepseek-v4-flash","headline":"This paper claims that tensor rational block Arnoldi and Lanczos algorithms reduce multidimensional linear time-invariant systems and solve large-scale Lyapunov tensor equations by projection onto rational Krylov subspaces.","keywords":["tensor rational Krylov subspace","Einstein product","model order reduction","multilinear dynamical systems","Lyapunov tensor equations","block Arnoldi algorithm","block Lanczos algorithm","adaptive interpolation points"],"falsifier":"Run Algorithm 4 on a small random stable tensor $A$ and random $B$ with $K_1=K_2=2$, form $X_m=\\mathcal{V}_m*Y_m*\\mathcal{V}_m^T$ from the reduced equation (6.5), compute the true residual $R_m=A*X_m+X_m*A^T+B*B^T$, and compare $\\|R_m\\|$ with the bound $2\\|\\Gamma_{m,A}*X_m*\\mathcal{V}_m^T\\|$ from Theorem 6.1; if the bound is violated, or if $\\mathcal{W}_m^T*(R_m*\\mathcal{W}_m)$ is not near zero, the Galerkin reduction is not valid as stated. Also check numerically whether $\\mathcal{V}_m*\\mathcal{W}_m^T*B$ equals $B$; for generic data it will not.","tokens_in":20255,"feed_emoji":"🧮","tokens_out":10324,"duration_ms":80384,"temperature":0.7,"pith_summary":"Multidimensional linear time-invariant systems, whose state, input, and output are tensors, grow too large for direct simulation and control design. This paper proposes to reduce such systems by projecting them onto tensor rational Krylov subspaces built with the Einstein product, and introduces two algorithms: tensor rational block Arnoldi and tensor rational block Lanczos. The claim is that these projections produce low-order tensor systems whose transfer functions approximate the original, with error bounds, and that the Lanczos variant also gives low-rank approximate solutions to the Lyapunov tensor equations needed for balanced truncation. An adaptive rule for choosing interpolation points from residual-based error estimates is part of the proposed method.","feed_headline":"Two tensor rational Krylov algorithms shrink multidimensional systems","feed_subtitle":"Projection with the Einstein product builds smaller models and solves Lyapunov tensor equations.","key_machinery":"The central object is the tensor rational block Krylov subspace $\\mathcal{K}_m(A,B,\\Sigma_m)=\\operatorname{Range}\\{(A-\\sigma_1 I)^{-1}*B,\\dots,(\\prod_{k=1}^m (A-\\sigma_k I)^{-1})*B\\}$, built one block at a time. In the Arnoldi variant each new block is $(A-\\sigma_j I)^{-1}*\\mathcal{V}_j$, orthogonalized against previous blocks by a tensor QR step; in the Lanczos variant the same operation is applied to $A^T$ with a second start tensor $C^T$, and bi-orthonormality is enforced by QR and SVD normalizations. The whole construction uses only the Einstein product and block-tensor notation, so the reduced system has the same tensor format as the original. These bases are what convert a large MLTI system or Lyapunov equation into a small projected problem.","core_discovery":"The central claim is that rational Krylov projection can be carried out entirely in tensor form, without unfolding the system into a matrix. For a system $\\dot{X}(t)=A*X(t)+B*U(t)$, $Y(t)=C*X(t)$, the paper builds bases $\\mathcal{V}_m$ and $\\mathcal{W}_m$ for coupled rational Krylov subspaces $\\mathcal{K}_m(A,B,\\Sigma)$ and $\\mathcal{K}_m(A^T,C^T,\\Sigma)$, forms reduced tensors $A_m=\\mathcal{W}_m^T*(A*\\mathcal{V}_m)$, $B_m=\\mathcal{W}_m^T*B$, $C_m=C*\\mathcal{V}_m$, and states that the reduced transfer function approximates the original to the accuracy controlled by the error bounds in Theorems 3.1–3.3. The same Lanczos bases are used to project the continuous Lyapunov tensor equation $A*X+X*A^T+B*B^T=O$ to a small equation $T_m*Y_m+Y_m*T_m^T+B_m*B_m^T=O$, with the residual norm bounded in Theorem 6.1 and a factored low-rank approximation proposed. The numerical examples with random sparse tensors and a discretized 2D heat equation are offered as evidence that the reduced transfer functions track the original frequency response.","pith_inferences":["A reader checking the proof of Theorem 6.1 will find the step $\\mathcal{V}_m*\\mathcal{W}_m^T*B=B$ used to eliminate $B*B^T$ from the residual; Algorithm 4 starts from $(A-\\sigma_1 I)^{-1}*B$ rather than $B$, so this identity is not automatic and the residual bound and the reduced Lyapunov equation may need an extra term.","The same tensor rational Krylov construction should extend to order $N>2$ by the block-tensor definitions in Section 2; the paper only runs experiments for 4th-order tensors but the algebraic relations are written for general $N$.","A natural next experiment is to isolate the benefit of the adaptive shift rule by comparing it to fixed shifts on the same heat-equation example, reporting transfer error versus total number of solves; the paper's tables compare the rational Lanczos method to the classical block Lanczos method but not to non-adaptive rational shifts.","The error estimates in Table 3.1 are heuristic simplifications of the exact residual expression; a practical implementation would need to test which of the six estimates tracks the true transfer-function error, since the paper's experiments do not report that comparison."],"forward_implications":["If the algorithms perform as claimed, a package that handles tensor arithmetic can reduce an MLTI system without ever forming the large unfolded matrix, saving memory.","The adaptive interpolation-point rule can be reused in both Arnoldi and Lanczos tensor settings, and it removes the need for a user-supplied set of shifts.","The Lyapunov solver produces a factored low-rank approximation $X_m \\approx Z_1*Z_2^T$, the form balanced truncation needs to compute Hankel singular values.","The error estimates in Theorems 3.1 and 3.3 give a stopping criterion: iterate until the residual estimate falls below a tolerance, as used in Algorithm 5."],"supporting_citations":[{"why":"Introduces the MLTI transfer function and Einstein-product system framework that the paper projects onto rational Krylov subspaces.","marker":"[9]"},{"why":"Defines multilinear time-invariant system theory and the tensor state-space form used throughout.","marker":"[7]"},{"why":"Provides the tensor block Arnoldi algorithm whose rational version is Algorithm 3.","marker":"[10]"},{"why":"Gives the adaptive rational block Lanczos-type algorithm for matrices that Algorithm 4 generalizes to tensors.","marker":"[4]"},{"why":"Supplies the adaptive rational block Arnoldi method for large-scale MIMO systems that motivates the tensor Arnoldi variant.","marker":"[1]"},{"why":"Defines the balanced-truncation setting and Lyapunov tensor equations the paper solves in Section 6.","marker":"[8]"},{"why":"Supplies the stable block Lanczos procedure that the tensor Lanczos biorthogonalization steps imitate.","marker":"[3]"}],"fun_headline_variants":["Tensor rational Krylov cuts multidimensional model order","Rational Krylov without matrix unfolding for MLTI systems","Tensor rational Krylov reduces multidimensional linear systems","Adaptive tensor Krylov for model order reduction of MLTI","Tensor Krylov projection for Lyapunov and system reduction"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing assumption in the Lyapunov part is that the input tensor $B$ lies in the low-dimensional range of the Lanczos basis $\\mathcal{V}_m$, so that $\\mathcal{V}_m*\\mathcal{W}_m^T*B=B$; Algorithm 4 begins with $(A-\\sigma_1 I)^{-1}*B$, so this is generally not true, and the reduced Lyapunov equation follows from the Galerkin condition only if that assumption is repaired.","fun_headline_variants_meta":{"raw":{"variants":["Tensor rational Krylov cuts multidimensional model order","Rational Krylov without matrix unfolding for MLTI systems","Tensor rational Krylov reduces multidimensional linear systems","Adaptive tensor Krylov for model order reduction of MLTI","Tensor Krylov projection for Lyapunov and system reduction"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000738,"raw_usage":{"total_tokens":3296,"prompt_tokens":943,"completion_tokens":2353,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":559,"completion_tokens_details":{"reasoning_tokens":2278}},"tokens_in":559,"tokens_out":2353,"duration_ms":15351,"temperature":1.0,"reasoning_tokens":2278,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-12T11:24:50.781454+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run Algorithm 4 on a small random stable tensor $A$ and random $B$ with $K_1=K_2=2$, form $X_m=\\mathcal{V}_m*Y_m*\\mathcal{V}_m^T$ from the reduced equation (6.5), compute the true residual $R_m=A*X_m+X_m*A^T+B*B^T$, and compare $\\|R_m\\|$ with the bound $2\\|\\Gamma_{m,A}*X_m*\\mathcal{V}_m^T\\|$ from Theorem 6.1; if the bound is violated, or if $\\mathcal{W}_m^T*(R_m*\\mathcal{W}_m)$ is not near zero, the Galerkin reduction is not valid as stated. Also check numerically whether $\\mathcal{V}_m*\\mathcal{W}_m^T*B$ equals $B$; for generic data it will not.","supporting_citations":[{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Introduces the MLTI transfer function and Einstein-product system framework that the paper projects onto rational Krylov subspaces."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Defines multilinear time-invariant system theory and the tensor state-space form used throughout."},{"cited_title":"El Guide, A","cited_arxiv_id":null,"evidence_quote":"Provides the tensor block Arnoldi algorithm whose rational version is Algorithm 3."},{"cited_title":"Barkouki, A","cited_arxiv_id":null,"evidence_quote":"Gives the adaptive rational block Lanczos-type algorithm for matrices that Algorithm 4 generalizes to tensors."},{"cited_title":"ABIDI, M","cited_arxiv_id":null,"evidence_quote":"Supplies the adaptive rational block Arnoldi method for large-scale MIMO systems that motivates the tensor Arnoldi variant."},{"cited_title":"Data-Driven Model Reduction for Multilinear Control Systems via Tensor Trains","cited_arxiv_id":"1912.03569","evidence_quote":"Defines the balanced-truncation setting and Lyapunov tensor equations the paper solves in Section 6."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies the stable block Lanczos procedure that the tensor Lanczos biorthogonalization steps imitate."}],"review_version":1}