{"id":"7ec87f4d-e85d-4aac-95a9-eda5c745bb14","arxiv_id":"2505.20379","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":7.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":4,"one_line_summary":"A softmax-plus-squared-rates reparameterization turns PH moment fitting into unconstrained optimization, achieving sub-percent fits for up to 20 moments with up to 100 phases.","lead":"This paper introduces a new parameterization of phase-type distributions that turns moment fitting into unconstrained optimization, enabling gradient descent to fit large distributions. It reports fitting up to 20 moments with up to 100 phases at sub-percent errors on a test set of phase-type distributions.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The reported success rates are computed only on moment targets sampled from PH distributions, so the abstract's claim about fitting real-world moment sequences is untested on infeasible targets and is the load-bearing gap.","rationale":"The re-parametrization in Section 4 is mathematically sound: Propositions 1-3 establish surjectivity onto the interior of the respective PH families, and the moment formulas are standard. The internal evidence (Section 7) supports the claim that, on PH-generated targets, the method frequently fits up to 20 moments within the stated tolerances. The weakest point is the gap between the evaluated targets and the intended application. All 1500 targets in Section 6.1.1 are PH-generated; thus every instance is feasible for some PH(n), and the optimization is essentially an inversion problem. For an arbitrary moment sequence from real data, feasibility is not guaranteed, and the method could either fail to converge to a good fit or return a PH whose moments are far from the targets, with no way to know that the target is infeasible. The paper's own conclusion (Section 10) lists only the minimal-size limitation, not this feasibility gap. Because the headline claim is an empirical generalization to 'given moments' / 'real-world data', the missing evaluation on infeasible or non-PH targets is load-bearing. The proposed test—applying the identical protocol to non-PH moment sequences—would settle whether the claim requires an explicit feasibility caveat. Since the reader's verdict already conditions on this point, the verdict should remain Conditional.","tokens_in":21474,"tokens_out":9000,"duration_ms":89883,"concrete_test":"Run the published implementation on a new test set of moment sequences that are not sampled from PH distributions: for example, the first 20 moments of lognormal, Weibull, and Pareto distributions (with first moment normalized to 1) or random Hankel-positive moment sequences with no known PH(n) representation. Use the same protocol as Section 6.2: n=100, l=20, η=0.5%, and the multi-start search with s=10000. Compare the success rate to Figures 5a-5c. If the non-PH success rate is close to the PH-generated rate, the concern is resolved; if it falls substantially (e.g., below 50%), the abstract's 'vast majority' claim should be explicitly restricted to PH-feasible targets.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central empirical claim (abstract; Section 7) rests entirely on the test set of Section 6.1.1, where every target moment sequence is generated by sampling a PH distribution (General, Coxian, or Hyper-Erlang) and computing its moments. For every such target a feasible PH representation is known to exist for sufficiently large n, so the optimization problem is guaranteed to have a solution with zero error in the large-n limit. The reported success rates therefore measure how often the heuristic search finds an approximate preimage on a feasible set. Real-world moment sequences, however, need not be PH-realizable; for such inputs there may be no PH(100) with all 20 moments within 0.5%, and the method has no feasibility check or certified error bound. Section 1 motivates the method precisely with real-world data, so the load-bearing premise of the headline claim is that the high success rates transfer to non-PH targets. This is not established, and the paper does not flag it as a limitation. The issue is not internal inconsistency; it is an external-validity gap in the central claim.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes an unconstrained optimization approach for fitting moment sequences with phase-type (PH) distributions. It introduces differentiable reparameterizations for general PH, Coxian, and Hyper-Erlang distributions, mapping unconstrained parameters onto valid Markovian representations. Moment matching is formulated as a weighted least-squares objective optimized by gradient descent with multiple random starts. The authors report experiments on 1500 moment targets sampled from PH distributions, claiming success in fitting up to 20 moments with PH distributions of up to 100 phases within relative errors below 0.5% in most cases. They also demonstrate joint fitting of moments and CDF quantiles, and a GI/GI/1 queueing case study.","tokens_in":21723,"tokens_out":11353,"duration_ms":112068,"significance":"If the central claims hold, the paper would be a useful practical advance: it scales moment fitting to more moments and larger PH sizes than existing exact or iterative methods, and it provides a simple, differentiable parameterization that can be reused for other differentiable fitting objectives. The propositions in Section 4 are conceptually important for the method, the GitHub implementation is a concrete strength, and the queueing application illustrates a plausible use case. However, the significance is tempered by two issues: the proof of the main general reparameterization proposition contains a technical error, and the empirical evaluation only uses targets sampled from PH distributions, so the headline claim about fitting 'given moments' from real-world data is not directly tested on non-PH or infeasible targets.","major_comments":[{"comment":"The converse direction of Proposition 1 is not proven correctly. The construction sets Z=log(E) for E=D+I+diag(d), but E has zero row sums and generally has negative diagonal entries, so log(E) is not defined and, even when it is, softmax(log(E)) is not equal to E. For example, a row of E of the form (2/3, 1/3, 1/3) with zero row sum would require softmax(log(E)) to equal E, which is false because softmax of a vector is a probability vector and E is not. A correct construction exists: define the stochastic matrix S by S_ij = T_ij/(-T_ii) for i ≠ j and S_ii = 1 - sum_{j ≠ i} S_ij, then set Z=log(S); this requires strict row-sum negativity, i.e., an interior condition that the proposition statement does not mention. Because Proposition 1 is the foundation of the general reparameterization, the proof must be corrected and the boundary condition stated precisely.","section":"Section 4.1, Proposition 1 and Eq. (5)"},{"comment":"The empirical evaluation only uses moment sequences generated by sampling PH distributions (General, Coxian, and Hyper-Erlang). For such targets a feasible PH representation is known to exist in the large-n limit, so the reported success rates measure how often the heuristic finds an approximate preimage on a feasible set. The abstract and Section 1 motivate the method with real-world data, but no experiment uses non-PH targets or infeasible moment sequences. As a result, the central empirical claim that the method 'is, in the vast majority of cases, capable of fitting' arbitrary prescribed moments is not established. I recommend adding experiments with targets sampled from non-PH families (e.g., lognormal, Weibull, gamma, or empirical moment estimates) and reporting success rates and the resulting errors, or explicitly restricting the claims to PH-realizable targets. The paper should also state whether any feasibility check is available for a user who only has a moment sequence.","section":"Section 6.1.1 and Section 7"},{"comment":"The experimental protocol reports the best result over multiple random starting points and, according to Section 5, over the fitted PH size n: 'we attempted to fit each sequence of moments in several trials with different hyperparameter settings, and extract and report the best fit. This includes the value n.' The success rates are therefore upper bounds relative to a single-run procedure, and it is unclear whether the bars in Figures 4-6 for each n are computed independently or as the best over n. This is load-bearing for the empirical claims, so the paper should clarify the protocol and, ideally, report the distribution of outcomes over starting points or at least the success rate for a single representative run.","section":"Section 5 and Section 6.2"}],"minor_comments":[{"comment":"Propositions 2 and 3 state that the images of the reparameterizations are the full sets of valid Coxian and Hyper-Erlang distributions, but the sigmoid and softmax mappings are strictly positive and therefore exclude boundary cases such as p_j = 0 and ω_j = 0. The statements should say 'interior' or 'arbitrarily close' in the same way as Proposition 1.","section":"Section 4.2 and Section 4.3"},{"comment":"The sampling procedure first draws the PH size uniformly from [1,200], but the fitting experiments use n ∈ {20,50,100}. It is not specified how targets with original size larger than the fitted size are treated, nor whether such targets are included in the reported success rates; this matters for interpreting results where the target may be infeasible for the fitted class.","section":"Section 6.1.1"},{"comment":"The objective uses weights w_i = m_i^{-2}, making it a sum of squared relative errors, while the success metric is the maximum relative error over all moments. The relation between the optimized objective and the reported success criterion is not discussed; a sentence explaining the choice would help.","section":"Section 6.2"},{"comment":"The KL-divergence values reported for the shape-fitting example are given without describing how the densities were discretized or how the KL divergence was computed; this should be stated for reproducibility.","section":"Section 8"},{"comment":"There are several minor typographical issues: 'preformed' should be 'performed' in Section 9, 'The rational' should be 'The rationale' in Section 8, and 'Kurtosios' in the caption of Figure 3b should be 'Kurtosis'.","section":"Section 9 and Figure 3b"}],"recommendation":"major_revision","confidential_remarks":"The paper presents a useful heuristic with a fixable proof issue and a clear external-validity gap. The proof of Proposition 1 can be repaired locally, and the addition of non-PH target experiments would substantially strengthen the empirical claims. I do not see an internal inconsistency that would warrant rejection, but the current abstract overstates what is demonstrated."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"You should look at this one. The core idea is simple and it works: reparametrize the PH representation as a = log(alpha), rates as gamma^2, and transition probabilities as row-softmax(Z), which turns the constrained moment-matching problem into unconstrained gradient descent. Propositions 1–3 are correct, the maps are differentiable, and they are surjective onto the interior of the relevant PH families. That is a real contribution, not a restatement of the existing exact-matching literature. Extending the same idea to Coxian and Hyper-Erlang with lower-dimensional parameter spaces is also sensible. The paper deserves credit for shipping a practical method that scales to 100 phases and 20 moments, which is beyond what the exact moment-matching methods can do.\n\nThe soft spot is the one the stress-test note flags, and it is real. Every target moment sequence in Section 6.1.1 is sampled from a PH distribution, so feasibility is guaranteed. The success rates therefore measure how often a nonconvex search finds a preimage on a feasible set. Real moment sequences need not be PH-realizable, and the paper never tests that case, nor does it provide a feasibility check or an error certificate. The abstract and conclusion phrase the claim as if it covers real-world data, which is a step too far. There are also smaller issues: the reported success rates are the best over multiple PH sizes and 10,000 starting points, which biases them upward; the code repository is mentioned but not pinned to a commit; and the queueing application is illustrative rather than a systematic study. None of this sinks the method. For a methods paper, the honest framing would be: here is a scalable search heuristic for PH moment fitting, it performs very well on feasible targets, and transfer to infeasible ones remains to be tested.\n\nThe citation pattern looks fair. The comparison to Buchholz and Kriege is brief but reasonable, and the appendix at least attempts a head-to-head. The shape-fitting section is a nice bonus, and the KL-divergence numbers in that example are plausible.\n\nWho gets value from this: stochastic modelers who currently stop at five moments because anything more is impractical, and researchers working on PH fitting generally. It deserves a serious referee. My recommendation: send it to review, but request that the claims be calibrated to the actual test set, and ideally add at least one experiment with moment sequences from non-PH distributions or explicitly discuss feasibility. That is a revision, not a rejection.","headline":"The unconstrained reparametrization is a clean and genuinely useful contribution, but the headline success rates are measured only on PH-generated targets, so real-world moment fitting is a plausible promise rather than a demonstrated fact.","tokens_in":22217,"tokens_out":903,"would_cite":true,"duration_ms":12524,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["90C26","60J27","62E17"],"pacs":[],"model":"deepseek-v4-flash","headline":"A re-parametrization of phase-type distributions turns moment fitting into unconstrained optimization, matching 20 moments with 100 phases to under 0.5% relative error.","keywords":["phase-type distributions","moment fitting","unconstrained optimization","gradient descent","reparametrization","Coxian distribution","Hyper-Erlang distribution","queueing analysis"],"falsifier":"Take a moment sequence that is known not to be PH-feasible, for example moments of a heavy-tailed distribution with a coefficient of variation so high that no PH(100) can match the first three moments, and run the method; if it still reports all moments within 0.5% relative error, the claim extends beyond PH-feasible targets, and if it fails or reports large errors, the paper's success guarantee is limited to feasible sequences.","tokens_in":21321,"feed_emoji":"📈","tokens_out":9629,"duration_ms":85990,"temperature":0.7,"pith_summary":"The paper claims that matching a phase-type (PH) distribution to a prescribed set of moments can be scaled well beyond prior limits by writing the constrained Markovian parameters (α, T) as differentiable functions of free variables, turning the fit into unconstrained optimization. On moment sequences sampled from PH distributions, the method matches up to 20 moments with PHs of up to 100 phases, keeping the worst relative error under 0.5% in the vast majority of instances. The same machinery also fits distribution shape through CDF points together with moments, and a queueing case study shows that using more fitted moments improves the steady-state queue length estimates of a GI/GI/1 queue. This matters because PH distributions are a standard way to bring non-exponential behavior into Markov-chain analysis, and moment matching is a practical route for obtaining them from data.","feed_headline":"New parametrization fits 20 moments to 100-phase PH distributions","feed_subtitle":"Matching many moments to phase-type models was limited to small sizes; this method lifts the barrier.","key_machinery":"The load-bearing object is the differentiable transformation π (with analogues φ and ψ for Coxian and Hyper-Erlang structures) that maps free parameters into the Markovian representation: α = softmax(a), T = diag(γ²)·[softmax(Z) − (I + softmax(Z)∘I)], where ∘ is elementwise multiplication. This map makes the feasible set of PH parameters the image of an unconstrained Euclidean space, so gradient descent on the weighted regression objective automatically stays feasible. The exponential (softmax/sigmoid) squashing functions guarantee non-negativity and the diagonal γ² factor controls the negative diagonal of T, while the softmax row normalization ensures non-positive row sums.","core_discovery":"The central discovery is a differentiable re-parametrization of PH distributions that removes all feasibility constraints. For a general PH(n) with initial vector α and subgenerator matrix T, the paper defines unconstrained parameters a ∈ R^n, γ ∈ (R\\{0})^n, Z ∈ $R^{{n×n}}$ and the mapping π given by α = softmax(a) and T = diag(γ²)[softmax(Z) − (I + softmax(Z)∘I)], whose image is exactly the set of valid PH(n) representations with no zero elements. Since π is differentiable, any differentiable objective, in particular the weighted squared relative moment error with weights w_i = $m_i^{{-2}}$, can be minimized by gradient descent in the unconstrained space, and every result maps back to a valid PH. The paper proves analogous re-parametrizations for Coxian and Hyper-Erlang subclasses, and reports numerical experiments in which the method fits 5, 10, and 20 moments to PHs of size 20, 50, and 100, with success rates near 100% in many settings and worst-case relative errors below 0.5% for the majority of targets. It also demonstrates joint fitting of moments and CDF quantiles, and a GI/GI/1 queueing study where the stationary queue length computed from PHs fitted to 2–5 moments converges to the ground truth as more moments are used.","pith_inferences":["Because the evaluation samples target moments only from PH distributions, the reported success rates do not yet cover moment sequences that are not PH-feasible; testing on empirical traces or on distributions outside the PH family would clarify how the method behaves when no exact PH(n) fit exists.","The re-parametrization covers only PH representations with no zero elements; boundary cases are approached but never reached, so problems whose optimal solution lies exactly on a boundary may need separate handling or a post-processing step.","The method's runtime (minutes to a few hours) fits offline model construction; a natural extension would be warm-starting or dimensionality reduction to make it usable in online or repeated-fitting settings, but the paper does not address this."],"forward_implications":["Moment fitting is no longer confined to small problems: PH distributions with dozens to hundreds of phases can be fitted to 20 moments, which was previously impractical.","Any differentiable property of a PH distribution (PDF, CDF, hazard rate, Laplace transform) can be added to the objective function, so moment fitting and shape fitting can be performed jointly in a single optimization.","For queueing models, the method provides a way to quantify the influence of the i-th moment on performance measures such as the steady-state queue length in a GI/GI/1 queue.","The Coxian and Hyper-Erlang re-parametrizations fit into the same framework, so practitioners can trade off parameter count against generality within one codebase."],"supporting_citations":[{"why":"Supplies the definition and the moment formula m_i = i!(−1)^i αT^{−i}1_n for PH distributions, on which the objective is built.","marker":"Horváth and Telek (2024)"},{"why":"Sets the prior state of the art for moment matching and the limits on the number of moments and PH size that the new method claims to overcome.","marker":"Johnson and Taaffe (1991, 1989)"},{"why":"The heuristic constrained-optimization approach to PH moment fitting that the paper compares against and whose scalability limitations motivate the unconstrained re-parametrization.","marker":"Buchholz and Kriege (2009)"},{"why":"Provides the four evaluation criteria (number of moments, efficiency, generality, minimality) and the closed-form Erlang–Coxian matching method that the new approach extends.","marker":"Osogami and Harchol-Balter (2006)"},{"why":"Demonstrates that moments of inter-arrival and service times suffice for predicting GI/GI/1 queues, motivating the practical need for high-order moment fitting.","marker":"Baron et al. (2024)"},{"why":"Supplies the Quasi-Birth-and-Death method used to compute the PH/PH/1 steady-state queue length in the application section.","marker":"Harchol-Balter, Mor (2013)"}],"fun_headline_variants":["Unconstrained optimization fits 20 moments to 100-phase PHs","Scale moment fitting: 20 moments, 100-phase PHs, no constraints","Reparametrize to unconstrained space, fit 20 moments and 100 phases","PH fitting without constraints: 20 moments up to 100 phases","Moment matching goes large: 20 moments, 100 phases, unconstrained"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The paper's headline accuracy numbers come from moment sequences that were themselves sampled from PH distributions, so a feasible PH representation is known to exist; for arbitrary real-world moment sequences, feasibility is not guaranteed and the method's success rate on such targets is untested.","fun_headline_variants_meta":{"raw":{"variants":["Unconstrained optimization fits 20 moments to 100-phase PHs","Scale moment fitting: 20 moments, 100-phase PHs, no constraints","Reparametrize to unconstrained space, fit 20 moments and 100 phases","PH fitting without constraints: 20 moments up to 100 phases","Moment matching goes large: 20 moments, 100 phases, unconstrained"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000423,"raw_usage":{"total_tokens":2268,"prompt_tokens":1140,"completion_tokens":1128,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":756,"completion_tokens_details":{"reasoning_tokens":1027}},"tokens_in":756,"tokens_out":1128,"duration_ms":10958,"temperature":1.0,"reasoning_tokens":1027,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-07T13:56:56.889252+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Take a moment sequence that is known not to be PH-feasible, for example moments of a heavy-tailed distribution with a coefficient of variation so high that no PH(100) can match the first three moments, and run the method; if it still reports all moments within 0.5% relative error, the claim extends beyond PH-feasible targets, and if it fails or reports large errors, the paper's success guarantee is limited to feasible sequences.","supporting_citations":[],"review_version":1}