{"id":"27d7060a-f58d-4b3a-a4ae-945c9d776192","arxiv_id":"2505.16906","paper_version":2,"verdict":"REJECT","confidence":"HIGH","novelty_score":2.0,"correctness_risk":"high","formal_verification":"none","parameter_count":2,"one_line_summary":"The higher order Jacobi method is a known repeated-squaring restatement of Jacobi iteration framed as a neural network, and the reported near-constant GPU scaling is not credible.","lead":"This paper rewrites the classic Jacobi iteration for solving Ax=b by squaring its update matrix, so that one 'higher order' step equals many ordinary steps, and then recasts the recurrence as a deep linear network. It also reports GPU timings that it reads as a dramatic complexity improvement, but the timing analysis is too fragile to support that claim.","discovery_kind":"incremental","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The claimed near-constant GPU scaling rests on fitted exponents from last two data points; reported n^-0.07 inverse is physically impossible and likely reflects kernel-launch overhead, not algorithmic complexity.","rationale":"The reader's weakest assumption is that wall-clock timings over a narrow size range, condensed into exponents from the last two points, represent asymptotic complexity. This is indeed the single most load-bearing concern: the paper's abstract and conclusion stake the entire contribution on a 'significant enhancement in the order of complexity' on GPU, and Table 2 is the only evidence offered. The fitted exponents violate elementary lower bounds (reading an n x n matrix costs Omega(n^2); a sparse matrix with O(n) nonzeros costs Omega(n)), so the measurements must be dominated by fixed overheads or truncated ranges. The ad hoc exclusion rule is disclosed in the text but makes the fitting procedure non-verifiable, and no error bars or repeated trials are reported. I therefore agree with the reader's rejection. Secondary mathematical issues, such as Eq. (18)'s incorrect stepping relation and the conflation of 2k with 2^k in Eq. (15), are genuine but less central: even if those were fixed, the empirical scaling claim would still be unsupported. Conversely, if the proposed large-scale benchmark reproduced exponents above 1.5 and confirmed HOJM is slower than 20 matmuls, the paper's own data would still not justify the claimed lower complexity class. The reader's REJECT verdict stands unchanged.","tokens_in":9653,"tokens_out":6082,"duration_ms":47384,"concrete_test":"Run a controlled benchmark on a single H100 (System 2): solve and invert for n = 512, 1024, 2048, 4096, 8192, 16384 with HOJM, using torch.cuda.Event to time only the 20 matrix-matrix products, and separately time one cuBLAS matmul of the same dense matrices. Fit the exponent over the largest five sizes after subtracting an empty-kernel latency measurement. If the exponent is not above 1.5, or if HOJM runs faster than 20x the single matmul, the reported sublinear scaling is an artifact.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central claim is the 'transformative' GPU complexity enhancement reported in Table 2 (Section 3). The exponents are fitted from the last two data points, with an ad hoc rule that drops the middle point of the last three if a sudden jump is observed. This procedure can manufacture arbitrarily small exponents. More fundamentally, the reported n^0.19 solve and n^-0.07 inverse exponents are inconsistent with any algorithm that must read the matrix: reading a dense n x n matrix is Omega(n^2), and reading a sparse matrix with nnz=O(n) is Omega(n). A negative exponent for inversion is impossible unless the measurements are dominated by fixed launch and overhead effects. No error bars or repeated runs are given, and the vertical range of plots is capped at 20 seconds, which can truncate the large-size behavior. Thus the load-bearing premise that HOJM achieves a lower order of complexity is unsupported; the data are more consistent with overhead-dominated timings over a narrow size range.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes a 'higher order Jacobi method' (HOJM) for solving linear systems and computing matrix inverses. The method recursively squares the Jacobi iteration matrix and accumulates a bias term, so that k higher-order steps are equivalent to 2^k ordinary Jacobi iterations. The resulting product formula is interpreted as a closed-form expression for the inverse of the normalized matrix, and the structure is likened to a deep linear network. The paper also proposes an 'induced source' iterative scheme for reusing a known inverse when the coefficient matrix changes slowly. The main empirical claim, based on GPU timing measurements, is that HOJM achieves a dramatically lower order of complexity than standard solvers, with fitted exponents as low as n^0.19 for solving and n^(-0.07) for inversion. The paper includes code on GitHub.","tokens_in":9792,"tokens_out":5196,"duration_ms":40745,"significance":"If the near-constant scaling claim were valid, it would indeed be a transformative result, as no algorithm that reads a matrix of size n can run in time sublinear in the matrix size. The paper's algebraic derivation of the higher-order iteration and the product formula for the inverse is elementary and is presented as a connection between iterative methods and deep linear networks. The availability of code and reproducibility scripts is a strength. However, the load-bearing empirical claim rests on exponent fitting from two data points, and the reported exponents are physically impossible for any algorithm that must read the input matrix. The paper's own Discussion section also acknowledges that HOJM has the same computational complexity as matrix-matrix multiplication, which contradicts the claimed 'significant enhancement in the order of complexity.' Without support for the central performance claim, the contribution reduces to an interesting but straightforward algebraic reformulation of Jacobi iteration.","major_comments":[{"comment":"Equation (18) states x(2k) = α(k)·x(2k−1) + β(k), claiming that the k-th order iterate can be obtained from the previous higher-order guess x(2k−1). This is not correct. Equation (15) gives x(2k) in terms of x(0) and the k-th order coefficients; applying the k-th order operator to x(2k−1) does not reproduce x(2k) because the k-th order operator advances by 2^k Jacobi steps, not by one step. The correct recurrence is x(m + 2^k) = α(k)·x(m) + β(k) for general m, or in the paper's indexing x(2k) = α(1)·x(2k−2) + β(1) for the second-order case. As written, Eq. (18) is mathematically inconsistent with Eq. (15) and would produce incorrect iterates if implemented. Since this equation is used to motivate the 'faster convergence' claim and the deep recurrent network analogy, it is a load-bearing error that must be corrected.","section":"2.1, Eq. (18)"},{"comment":"The empirical complexity exponents are computed using only the last two data points, with an ad hoc rule that drops the middle of the last three points if a 'sudden jump' is observed. As stated in Section 3: 'The runtime empirical computational complexity for each curve is calculated using the last two data points. If a sudden jump is observed in the last 3 data points, the middle point is excluded for the calculation.' This procedure can produce arbitrarily small exponents. The reported exponents n^0.19 for solving and n^(-0.07) for inversion are below the O(n^2) cost of even reading a dense n×n matrix, and a negative exponent for matrix inversion is impossible for any algorithm that must output an n×n result. These exponents therefore reflect fixed kernel-launch overhead and memory-bound effects over a narrow size range, not algorithmic complexity. The absence of error bars, repeated runs, or any discussion of measurement variability further undermines the claim. This is the central load-bearing evidence for the 'significant enhancement in the order of complexity' stated in the abstract.","section":"Section 3, Table 2"},{"comment":"The Discussion states that 'The proposed higher-order method has the same computational complexity as matrix-matrix multiplication' (naively O(n^3) for dense matrices and O(n^3) for the product of n×n matrices). This directly contradicts the abstract's claim of 'significant enhancement in the order of complexity.' The GPU timings in Table 2 indicate only that, over the measured size range, the wall-clock time grows slowly; they do not establish a lower order of complexity. The paper conflates measured wall-clock time on a specific GPU implementation with asymptotic algorithmic complexity. At minimum, the paper must clearly separate the algorithmic complexity of HOJM from the observed speedup due to GPU parallelization and must temper the 'order of complexity' language accordingly.","section":"Section 4 (Discussion) vs. Abstract"},{"comment":"The induced source iterative scheme is asserted to converge for slowly varying coefficient matrices, but no convergence proof or numerical demonstration is provided. The scheme involves the free relaxation parameter ω, and the iteration (27)–(30) is a fixed-point iteration whose convergence depends on properties of the unknown function 𭟋 and the matrix difference A − A_d. No conditions on ω, on the rate of variation of A_d, or on the spectral radius of the iteration operator are given. Since this section is the basis for the claimed ability to 'efficiently resolve system variations without recomputing the coefficients,' the absence of any analysis or experiment makes this claim unsupported.","section":"2.4, Eqs. (27)–(30)"}],"minor_comments":[{"comment":"The paper uses two different acronyms for the proposed method: HOJM in the abstract and Section 1, and 'HOM' in the second paragraph of Section 1. Use one consistently.","section":"Throughout"},{"comment":"The product notation 0∏i=k−1 is nonstandard and the order of matrix multiplication matters. It should be written as ∏_{i=k−1}^{0} (α(i) + I) with the convention that the product is taken in reverse order, or as an explicit descending product.","section":"2.2, Eq. (22)"},{"comment":"The table contains two rows labeled 'after' (after System 1 and after System 2) whose meaning is unclear. These rows appear to belong to the HOJM columns but are not explained in the text. Clarify or remove.","section":"Table 2"},{"comment":"There are typographical errors: 'single presision' should be 'single precision,' and 'coeffecient' should be 'coefficient.' Also, 'The best order in each row has be highlighted' should be 'has been highlighted.'","section":"Section 3"},{"comment":"The text states that the vertical range is capped to 20 seconds, but for log-log plots the cap is a horizontal truncation of large-size behavior. The paper should justify this cap and show whether the truncated points would change the fitted exponents.","section":"Figure 1"}],"recommendation":"reject","confidential_remarks":"The paper is posted in cond-mat.supr-con, but its content is a numerical linear algebra and GPU computing study; this is a scope mismatch for the journal. The GitHub code and reproducibility effort are commendable, but the central empirical claim is not supportable. The reported exponents (n^0.19 and n^(-0.07)) violate information-theoretic lower bounds and are clearly an artifact of the last-two-points fitting rule. The algebraic error in Eq. (18) is also a serious flaw. Even with a corrected Eq. (18), the paper would need to be rewritten to present the method as an alternative GPU-accelerated iterative technique rather than as an algorithm with a fundamentally lower order of complexity. I see no incremental revision that would salvage the current claims within the manuscript's scope."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Colleague,\n\nThe short version: the central claim of a transformative GPU speedup does not survive contact with the data. The reported exponents, n^0.19 for solving and n^-0.07 for inverse, are below the cost of reading the matrix, so they almost certainly reflect fixed kernel-launch overhead across a narrow size range, not algorithmic complexity. The algebraic core is real but classical: the recursion α(k)=α(k-1)^2 and β(k)=(α(k-1)+I)β(k-1) is the repeated-squaring identity for (I-α)^{-1}=∏(I+α^{2^i}), which requires ρ(α)<1.\n\nWhat is actually new? The induced source scheme in Section 2.4, which tries to reuse a known inverse when the coefficient matrix varies slowly, is a genuine proposal. But it comes with no convergence analysis and no numerical tests, so it is a suggestion, not a result. The neural-network framing is a relabeling: a linear recurrence is not a deep network in any meaningful sense.\n\nWhat the paper does well: the derivation up to Eq. (17) is algebraically correct, and the authors are honest that the method inherits the convergence conditions of classical Jacobi. They also ship code and data on GitHub, which is more than many submissions.\n\nThe soft spots are serious. Eq. (18) states x(2k)=α(k)x(2k−1)+β(k), but the correct predecessor is x(2k−2). Eq. (15) conflates 2k and 2^k. The empirical complexity in Table 2 is fitted from the last two data points, with an ad hoc rule that drops the middle of the last three if there is a jump; no error bars or repeated runs are reported. The negative exponent for inversion is impossible for any algorithm that touches the input. The 20-second cap on the y-axis truncates large-size behavior, which can hide where the asymptotics actually kick in.\n\nThese flaws are not minor. They undermine the paper's headline claim. But the underlying identity is correct, and the induced-source idea might be worth a paragraph somewhere. As it stands, I would not cite this paper, and I would not send it to a serious journal. A desk reject is defensible. If it ever comes back with real benchmarks, a fixed Eq. (18), and a convergence proof for the induced-source scheme, it could be a small but useful note.\n\nFor peer review: I lean toward accept-to-review only if the editor believes the induced-source idea has merit; otherwise, reject.","headline":"Repeated-squaring Jacobi is a known identity; the GPU scaling exponents are physically implausible and the empirical methodology is ad hoc, but the induced-source proposal is worth a footnote.","tokens_in":10377,"tokens_out":3931,"would_cite":false,"duration_ms":29381,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":false},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["65F10","65F05","15A09"],"pacs":[],"model":"deepseek-v4-flash","headline":"This paper claims that recursively squaring the Jacobi update produces a finite product identity for the inverse of the normalized system matrix, turning linear solves into GPU-friendly matrix chains.","keywords":["higher order Jacobi method","matrix inverse","iterative linear solvers","matrix-matrix multiplication","GPU acceleration","deep linear networks","physics-informed neural networks","computational complexity"],"falsifier":"Measure HOJM-GPU runtime over at least a factor of 16 in matrix size, fit the log-log exponent over all points, and compare with the cost of the matrix-matrix multiplications the method performs; the near-constant scaling claim is false if the fitted exponent rises toward the matrix-multiplication exponent rather than staying near zero.","tokens_in":9384,"feed_emoji":"🧮","tokens_out":12746,"duration_ms":96371,"temperature":0.7,"pith_summary":"The paper claims that the classical Jacobi iteration can be reorganized into a higher-order method: by recursively squaring the iteration matrix, one step at order $k$ advances the solution as far as $2^k$ ordinary Jacobi steps, and the same recurrence yields an explicit finite product for the inverse of the normalized coefficient matrix, $A_n^{-1} = \\prod_{i=0}^{k-1}(\\alpha^{(i)}+I)$. If that identity holds, solving a linear system becomes a chain of matrix-matrix multiplications rather than many matrix-vector sweeps, which is exactly the operation GPUs execute fastest. The paper further proposes an induced-source iteration that reuses a precomputed inverse to solve nearby systems with only matrix-vector work, and it presents the whole construction as a training-free deep linear network with analytically defined weights. On two GPU systems, the authors report runtime scaling exponents for the higher-order method far below those of built-in dense solvers, and they take this as evidence that the method is significantly faster in practice.","feed_headline":"Recursive squaring turns Jacobi sweeps into one inverse product","feed_subtitle":"A finite product of squared iteration matrices computes the inverse directly and scales almost flatly on GPUs.","key_machinery":"The central object is the recursively squared iteration matrix $\\alpha^{(k)} = (\\alpha^{(k-1)})^2$ with the companion bias $\\beta^{(k)} = (\\alpha^{(k-1)} + I)\\beta^{(k-1)}$, seeded by $\\alpha^{(0)} = -(A_n - I)$ and $\\beta^{(0)} = b_n$. The carrying identity is Eq. (25), $A_n^{-1} = \\prod_{i=0}^{k-1}(\\alpha^{(i)} + I)$, which packages $2^k$ Jacobi sweeps into one matrix-product chain and is what turns the iterative method into a direct inverse formula. The induced-source update, Eqs. (27)-(30), is the companion mechanism that reuses the precomputed inverse by iterating only the source term when the coefficient matrix changes slowly.","core_discovery":"On its own terms, the central discovery is that the Jacobi recurrence $x^{(k+1)} = \\alpha x^{(k)} + \\beta$ can be exponentiated: defining $\\alpha^{(k)} = (\\alpha^{(k-1)})^2$ and $\\beta^{(k)} = (\\alpha^{(k-1)} + I)\\beta^{(k-1)}$ makes one order-$k$ update equivalent to $2^k$ classical sweeps, and the expanded bias term gives the closed form $A_n^{-1} = \\prod_{i=0}^{k-1}(\\alpha^{(i)} + I)$ for the inverse of the normalized matrix. The paper treats this equation as an explicit finite construction of the inverse, so that direct and iterative methods are unified rather than opposed. It also reads the coefficient hierarchy as a deep linear network with weights and biases fixed by the physics, bypassing the training step of a physics-informed neural network. For slowly changing systems, the proposed induced-source scheme (Eqs. (27)-(30)) keeps the precomputed inverse and iterates only on the source vector. The empirical section then reports GPU timings whose fitted exponents are $n^{0.19}$/$n^{0.79}$ for solving and $n^{0.24}$/$n^{-0.07}$ for inversion, which the authors attribute to the parallel nature of matrix-matrix multiplication.","pith_inferences":["Beyond the paper: the same recursive squaring idea should apply to any stationary iteration whose iteration matrix is contractive, yielding analogous product formulas for Gauss-Seidel or successive over-relaxation.","Beyond the paper: the induced-source scheme could be combined with time stepping or Newton linearization so one precomputed inverse serves many successive linear systems in nonlinear or time-dependent problems, an extension the paper only gestures at.","Beyond the paper: because the higher-order products densify as they are squared, a concrete test would be to truncate each $\\alpha^{(i)}$ sparsely and measure how the residual of Eq. (25) degrades, quantifying the pruning-versus-accuracy trade-off in the same terms as runtime."],"forward_implications":["One application of an order-$k$ coefficient pair advances the solution as far as $2^k$ ordinary Jacobi steps, so the paper's choice of $k=20$ corresponds to about one million classical iterations.","The product identity $A_n^{-1} = \\prod_{i=0}^{k-1}(\\alpha^{(i)} + I)$ gives an explicit finite expression for the inverse of the normalized system matrix, so the method can be used as a direct solver whenever the Jacobi convergence conditions hold.","The recurrence defines all network weights $\\alpha^{(i)}$ and biases $\\beta^{(i)}$ analytically, meaning the linear solve needs no training phase, unlike a physics-informed neural network.","With a slowly varying coefficient matrix, the induced-source scheme (Eqs. (27)-(30)) performs only matrix-vector multiplications after the initial inverse is computed, so re-solving a nearby system avoids recomputing higher-order coefficients.","The reported GPU runtime exponents for the higher-order method are $n^{0.19}$ and $n^{0.79}$ for solving and $n^{0.24}$ and $n^{-0.07}$ for inversion, which the authors attribute to the parallelism of matrix-matrix multiplication."],"supporting_citations":[{"why":"Supplies the classical Jacobi iteration that the higher-order method reorganizes.","marker":"[18]"},{"why":"Provides the Jacobi-type iteration for matrix equations that motivates the generalization.","marker":"[19]"},{"why":"Gives the diagonal-dominance convergence condition that the closed-form inverse relies on.","marker":"[28]"},{"why":"Supplies a verified convergence result for Jacobi, grounding the equivalence claim's necessary conditions.","marker":"[30]"},{"why":"Defines the physics-informed neural network framework that the paper contrasts with its explicit weights.","marker":"[23]"},{"why":"Provides the Strassen matrix-multiplication baseline used to frame the method's complexity.","marker":"[41]"},{"why":"Gives the Coppersmith-Winograd complexity baseline for matrix multiplication.","marker":"[43]"},{"why":"Lists multigrid's linear complexity as the standard to which the method compares.","marker":"[40]"}],"fun_headline_variants":["Higher-order Jacobi: inverse via matrix products, no training","Jacobi iteration squared: one product equals exponential sweeps","Neural-network-inspired Jacobi computes inverse directly","Exponentiating Jacobi: fast inverse on GPUs, no training","Squaring Jacobi sweeps: direct inverse from matrix-matrix products"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The empirical result stands on taking the slope of GPU wall-clock time over a small set of matrix sizes as the asymptotic complexity of the method, even though the reported exponents are far below the cost of reading the matrix once.","fun_headline_variants_meta":{"raw":{"variants":["Higher-order Jacobi: inverse via matrix products, no training","Jacobi iteration squared: one product equals exponential sweeps","Neural-network-inspired Jacobi computes inverse directly","Exponentiating Jacobi: fast inverse on GPUs, no training","Squaring Jacobi sweeps: direct inverse from matrix-matrix products"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000392,"raw_usage":{"total_tokens":2091,"prompt_tokens":1004,"completion_tokens":1087,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":620,"completion_tokens_details":{"reasoning_tokens":1000}},"tokens_in":620,"tokens_out":1087,"duration_ms":8717,"temperature":1.0,"reasoning_tokens":1000,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-07T14:53:15.925168+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Measure HOJM-GPU runtime over at least a factor of 16 in matrix size, fit the log-log exponent over all points, and compare with the cost of the matrix-matrix multiplications the method performs; the near-constant scaling claim is false if the fitted exponent rises toward the matrix-multiplication exponent rather than staying near zero.","supporting_citations":[{"cited_title":"Iterative methods for solving matrix equations","cited_arxiv_id":null,"evidence_quote":"Supplies the classical Jacobi iteration that the higher-order method reorganizes."},{"cited_title":"The jacobi and gauss–seidel-type iteration methods for the matrix equation axb= c","cited_arxiv_id":null,"evidence_quote":"Provides the Jacobi-type iteration for matrix equations that motivates the generalization."},{"cited_title":"Convergence of matrix iterations subject to diagonal dominance","cited_arxiv_id":null,"evidence_quote":"Gives the diagonal-dominance convergence condition that the closed-form inverse relies on."},{"cited_title":"Verified cor- rectness, accuracy, and convergence of a stationary iterative linear solver: Jacobi method","cited_arxiv_id":null,"evidence_quote":"Supplies a verified convergence result for Jacobi, grounding the equivalence claim's necessary conditions."},{"cited_title":"Im- plementation of strassen’s algorithm for matrix multiplication","cited_arxiv_id":null,"evidence_quote":"Provides the Strassen matrix-multiplication baseline used to frame the method's complexity."},{"cited_title":"On the asymptotic complexity of matrix multiplication","cited_arxiv_id":null,"evidence_quote":"Gives the Coppersmith-Winograd complexity baseline for matrix multiplication."},{"cited_title":"Multigrid methods","cited_arxiv_id":null,"evidence_quote":"Lists multigrid's linear complexity as the standard to which the method compares."}],"review_version":1}