{"id":"1424282f-47c9-4467-9c19-5054c1f81988","arxiv_id":"2510.02664","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":5.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":2,"one_line_summary":"HOMC is a MATLAB package for higher-order Markov chains that computes transition power tensors, limiting distributions, ever-reaching probabilities, and mean first passage times.","lead":"This paper presents HOMC, a free MATLAB package for computing transition probabilities, long-run behavior, and first passage times in higher-order Markov chains. It is the first dedicated toolbox for these computations, though the mathematics it implements comes from the author's earlier papers.","discovery_kind":"extension","skeptic_critique":{"model":"deepseek-v4-flash","headline":"erp stopping rule has no tail error bound; slowly mixing chains can be misclassified with default tolerance","rationale":"I read the paper as a software presentation. The central claim is that HOMC can compute all important quantities reliably. The bpow, rcmat, mfptd, and mfpti routines are backed by definitions or theorems in [8,9,11]; the only new algorithmic decision is erp's stopping rule. The reader's weakest assumption identifies this exactly. I considered whether the 'first of its kind' claim is a problem, but that is a historical novelty claim not a correctness issue. The lack of tests/license is a reproducibility concern but secondary. The erp issue can cause concrete misclassification errors; it is therefore the most load-bearing. My suggested test would settle whether the concern actually lands; if erp's default output matches the exact solution to tolerance in slowly mixing chains, the concern is mitigated, but as written the paper provides no such evidence.","tokens_in":11683,"tokens_out":4923,"duration_ms":56056,"concrete_test":"Take a second-order chain (n=2 or 3) whose reduced first-order matrix Q has a sub-dominant eigenvalue 1-ε with ε=1e-4, e.g., by perturbing a regular chain toward reducibility. Compute erp(P,1e-6) and erp(P,1e-14), and also compute the exact ever-reaching tensor by making each target state absorbing in Q and solving the linear system (I-Q_sub)x=b for absorption probabilities. If the default-tol erp output differs from the exact values by more than 1e-3 in any entry, the stopping rule is unsafe for classification. This can be automated with a short MATLAB script.","verdict_should_be":"UNCHANGED","load_bearing_attack":"Section 3.1 defines the ever-reaching probability tensor F = Σ_{k=1}^∞ F^{[k]} (Eq. 5) and the erp function stops when max(abs(F^{[k]})) < tol, returning the partial sum. Since all entries are nonnegative, this is max(F^{[k]}) < tol. The omitted tail T_k = Σ_{l>k} F^{[l]} is not controlled. For finite-state chains, the tail after the first time a term falls below tol can be as large as tol/(1-ρ), where ρ is the spectral radius of the substochastic transition matrix for paths avoiding the target state. For nearly decomposable chains (ρ close to 1, e.g., 0.9999), the tail can be ~1e-2, far above tol. This directly threatens the state classification routine described in the same section, where a state is recurrent iff f_{iii3...im}=1 for all histories. A computed value 0.99 instead of 1 would misclassify a recurrent state as transient, or a transient state with true probability 0.9995 might be rounded to 1. No error bound is stated, and the default tol=1e-6 gives no assurance. This is the most fragile load-bearing premise because every other routine either reduces to standard linear algebra or has a convergence theorem cited from [9]; erp's heuristic cut-off is new to this software.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The manuscript presents HOMC, a MATLAB package for higher-order Markov chains. Its functions implement a tensor 'box' product and powers (bprod, bpow), linear indexing (lind), diagonal tensors (diagt), matricization/tensorization (t2mat, mat2t), construction of the reduced first-order chain (rcmat, krprod), ever-reaching probabilities (erp), and mean first passage times (mfptd, mfpti). The paper claims these tools can compute the k-step transition tensor, limiting distributions, ever-reaching probability tensors, and MFPT tensors, and can check ergodicity/regularity and classify states as recurrent or transient. The manuscript includes worked examples, including a residual check of the MFPT equation and first-order reductions to standard matrix results.","tokens_in":11962,"tokens_out":5122,"duration_ms":50795,"significance":"If the package performs as claimed, it would fill a genuine gap: no other dedicated MATLAB toolkit for higher-order Markov chains is known to the community, and the tensor-based formulation is natural for multidimensional transition data. The paper's internal checks are a real strength: Example 3.2 verifies the MFPT equation to zero residual, and Examples 2.4 and 3.3 reproduce standard first-order matrix results. The package builds on the authors' prior theoretical work [8,9,11], which is appropriate for a software presentation. The main weakness is the erp truncation rule, which has no tail error bound and directly affects the state-classification feature; additionally, the 'first of its kind' claim is not benchmarked against existing software, and the manuscript does not include a test suite or code listing.","major_comments":[{"comment":"The stopping test max(abs(F[k])) < tol does not control the omitted tail sum_{l>k} F[l]. Since all entries are nonnegative, a small per-step increment does not imply a small remaining tail; for a nearly decomposable chain the tail can be as large as tol/(1−ρ), where ρ is the spectral radius of the substochastic matrix for transitions avoiding the target state. With ρ close to 1, the tail can exceed tol by orders of magnitude. This is load-bearing because the same section uses the exact condition f_{ii3...im}=1 to classify recurrent states. Please provide a tail bound, or compute F by solving the linear tensor equation F = P + (F−F_d)⊠P (analogous to Eq. (7)), or at least clearly document the limitation and require the user to supply a much stricter tol after an explicit convergence check.","section":"§3.1, Eq. (5) and erp stopping rule"},{"comment":"The recurrent/transient classifier hinges on exact equality f_{ii3...im}=1. With the approximate erp output, a recurrent state may be returned as 0.999997 and a transient state with true probability 0.99999 may be rounded to 1. The default tol=1e-6 gives no assurance in either direction. This is not a purely theoretical concern: the classification routine is a headline feature of the package. The manuscript should either return the approximate probability together with a confidence statement, or use an exact method for the binary recurrent/transient decision.","section":"§3.1, state classification"},{"comment":"The paper's verification consists of four examples, but no automated test suite, code listing, or systematic function-reference mapping is included. For a software paper, this is insufficient to establish that the package behaves correctly beyond the displayed inputs. I recommend adding a small regression test suite and a table listing each function, its signature, and the examples that exercise it. This would also make the claim that HOMC computes 'all important quantities' in the cited works independently checkable.","section":"General, §2–§3, reproducibility"}],"minor_comments":[{"comment":"'First of its kind' is a strong claim. Either report a brief search of existing MATLAB toolboxes for Markov chains/HMMs/tensor methods, or soften the wording to 'to our knowledge, no dedicated package...'.","section":"Abstract and §1"},{"comment":"When defining I, specify explicitly that it is a left identity for the box product but not a right identity. This prevents confusion with standard tensor identity tensors and explains the non-associativity warning.","section":"§2.2, identity tensor"},{"comment":"The text states that erp used the first 67 terms, but the exact MATLAB call and formatting are not shown. Include the full command sequence and the MATLAB version/rounding behavior so that the displayed outputs can be reproduced bit-for-bit.","section":"Example 3.1"},{"comment":"The iterative stopping criterion is standard, but no error estimate is given for the fixed-point iteration. Please state whether the convergence proof in [9] provides a contraction rate that would justify interpreting tol as an a posteriori bound on the error.","section":"§3.2, mfpti"}],"recommendation":"major_revision","confidential_remarks":"The erp tail issue is the most important technical concern. I would not be comfortable accepting a software paper whose state-classification feature can silently misclassify recurrent states as transient or vice versa. The fix is local—add a tail bound or solve the exact linear equation—so I see this as a major-revision point rather than a rejection. Please also verify the 'first of its kind' claim against existing software before publication."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Short version: the package is real, the examples check out, and the paper is honest about what it implements. The one thing that needs fixing before I would rely on it is the erp stopping rule.\n\nWhat's actually new: a dedicated MATLAB toolbox for higher-order Markov chains, built around a box product that is different from other tensor products. The functions cover the expected quantities — k-step transitions, limiting distributions, reduced first-order chain, MFPT, erp, state classification. The implementation uses linear indexing with ind2sub to avoid materializing index arrays, which is a sensible memory choice. The examples are internally consistent: Example 3.2's residual is zero, and the first-order reductions match standard matrix results. So the package appears to do what it says for the cases shown.\n\nSoft spots: (1) The erp function sums terms until the latest term falls below tol; there is no bound on the omitted tail. For a nearly decomposable chain the tail after that point can be orders of magnitude larger than tol, which means the returned F and any state classification built on it can be wrong. The stress-test arithmetic is right. This is the one place where a new numerical method is introduced without analysis. (2) The package is not self-contained in the paper: no test suite, no version hash, no license, no mapping from printed outputs to the commands that produced them. A referee cannot verify the code from the text alone. (3) The 'first of its kind' claim is plausible but not checked against other tensor toolkits; that is minor.\n\nNone of these change the verdict. The math is cited from prior papers; the new content is the software. The package is worth having for a specialized community, and the paper is worth refereeing — but the referee should be instructed to look at the erp tail question and the code packaging. If those two are addressed, this is a useful contribution.","headline":"A legitimate software paper with a real artifact, but the erp stopping rule needs an error bound before the state-classification output can be trusted.","tokens_in":12440,"tokens_out":1952,"would_cite":false,"duration_ms":38747,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["60J10","15A69"],"pacs":[],"model":"deepseek-v4-flash","headline":"The paper introduces HOMC, the first dedicated MATLAB package for higher-order Markov chains, and shows how to compute their key tensors and classifications.","keywords":["higher order Markov chains","MATLAB","tensors","transition probabilities","ever-reaching probabilities","limiting probability distributions","mean first passage times","classification of states"],"falsifier":"Run erp on the two-state first-order chain with transition matrix P = [[1−α, 0.5], [α, 0.5]] for α = 10^−10 and default tolerance 10^−6. The routine stops when max(|F[k]|) < 10^−6, which for this chain happens around k = 21; the returned F_{12} is about 2×10^−9, while the true ever-reaching probability from state 2 to state 1 is 1 (the sum of the long-tailed geometric series). This divergence would show the truncation criterion is not a valid error bound.","tokens_in":11535,"feed_emoji":"🧮","tokens_out":10342,"duration_ms":80070,"temperature":0.7,"pith_summary":"The paper introduces HOMC, a MATLAB package for higher-order Markov chains—chains whose next state depends on several past states, not just the present one. Its claim is that this is the first dedicated software of its kind, and that it can compute the main numerical quantities from the recent theory of such chains: k-step transition tensors, limiting probability distributions, ever-reaching probability tensors, and mean first passage time tensors, plus regularity and ergodicity checks and state classification. If that is right, an entire toolkit of tensor-based Markov chain results becomes directly usable for numerical experiments and applications.","feed_headline":"First MATLAB package automates higher-order Markov chain analysis","feed_subtitle":"Computes transition tensors, limiting laws, hitting probabilities, and first-passage times from a single tensor input.","key_machinery":"The box product ⊠ of two m-th order n-dimensional tensors A and B is the tensor C with c_{i1...im} = Σ_j a_{i1 j i2 ... i_{m-1}} b_{j i2 ... im}. It reduces to ordinary matrix multiplication when m=2, but is non-associative for m≥3, which is why the package provides dedicated powering (bpow) rather than relying on matrix-style identities. This product is the algebraic engine: tensor powers give k-step transition probabilities, and the recurrences for first-passage and mean-first-passage tensors are box-product equations. Memory efficiency comes from generating index tuples on the fly with MATLAB's ind2sub rather than storing all index combinations.","core_discovery":"The central claim is that a single MATLAB package, built around a non-associative tensor product called the box product (⊠), can replace bespoke code for almost all numerical work with higher-order Markov chains. For an m-th order, n-dimensional transition tensor P, the box product defines tensor powers whose entries are exactly k-step transition probabilities; the same operation yields the ever-reaching probability tensor as a series of box products and the mean first passage time tensor as the solution of a tensor linear equation. The package also constructs the reduced first-order chain's transition matrix via a Khatri-Rao product, so limiting distributions can be obtained from a dominant","pith_inferences":["The truncation rule in the ever-reaching probability routine is a practical heuristic, not a certified numerical method; for chains with heavy-tailed first-passage times, the returned tensor can be wrong by far more than the tolerance. A tail bound or an alternative linear-system solve would make this routine reliable.","The non-associativity of the box product is a genuine conceptual difference from matrix Markov chains: standard matrix-power shortcuts do not transfer, so users must think in terms of the package's bpow rather than ordinary algebra.","The dense O(n^m) storage limits practical state spaces; extending the package to exploit sparse or structured transition tensors would open up higher-order chain models with much larger state spaces.","The routines assume the theory in the related papers (regularity/ergodicity for uniqueness); a useful diagnostic addition would be warnings when a chain fails the ergodicity check before computing mean first passage times."],"forward_implications":["Users can compute the limiting distribution of a higher-order chain by forming the reduced first-order transition matrix and taking a dominant eigenvector, without manual bookkeeping.","The ever-reaching probability tensor gives a direct state classification: a state is recurrent iff its diagonal entries in F are all 1, transient otherwise.","The package makes the recent theoretical results on higher-order chains available as off-the-shelf routines, lowering the barrier to numerical experimentation.","Because it also handles first-order chains (m=2) as a special case, it can serve as a single tool for a broad range of Markov chain models.","The O(n^m) memory footprint is the theoretical minimum for dense tensors, so the package scales as far as the underlying data structure can reasonably go."],"fun_headline_variants":["MATLAB package tames higher-order Markov chains","First MATLAB tool for higher-order Markov analysis","Compute Markov chain tensors with one MATLAB package","Box product powers higher-order Markov chain math","Higher-order Markov chains: now in MATLAB"],"cache_read_input_tokens":2304,"weakest_assumption_plain":"The ever-reaching probability function stops summing the infinite series as soon as the largest term falls below a tolerance, without a proof that the omitted tail is that small; for chains with slowly decaying first-passage probabilities this can silently misreport the ever-reaching probability by far more than the tolerance.","fun_headline_variants_meta":{"raw":{"variants":["MATLAB package tames higher-order Markov chains","First MATLAB tool for higher-order Markov analysis","Compute Markov chain tensors with one MATLAB package","Box product powers higher-order Markov chain math","Higher-order Markov chains: now in MATLAB"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.00077,"raw_usage":{"total_tokens":3203,"prompt_tokens":657,"completion_tokens":2546,"prompt_tokens_details":{"cached_tokens":256},"prompt_cache_hit_tokens":256,"prompt_cache_miss_tokens":401,"completion_tokens_details":{"reasoning_tokens":2492}},"tokens_in":401,"tokens_out":2546,"duration_ms":15094,"temperature":1.0,"reasoning_tokens":2492,"cache_read_input_tokens":256,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-04T12:38:43.920573+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run erp on the two-state first-order chain with transition matrix P = [[1−α, 0.5], [α, 0.5]] for α = 10^−10 and default tolerance 10^−6. The routine stops when max(|F[k]|) < 10^−6, which for this chain happens around k = 21; the returned F_{12} is about 2×10^−9, while the true ever-reaching probability from state 2 to state 1 is 1 (the sum of the long-tailed geometric series). This divergence would show the truncation criterion is not a valid error bound.","supporting_citations":[],"review_version":1}