REVIEW 3 major objections 5 minor 22 references
An Unconstrained Optimization Approach to Moment Fitting with Phase Type Distributions
T0 review · 3 major / 5 minor · reviewed 2026-08-07 · deepseek-v4-flash
Pith's one-line read 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.
desk verdict 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. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
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.
What would settle it
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.
Extended reading notes
Core claim
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.
Load-bearing premise
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.
Editorial extensions
If this is right
- 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.
Reading between the lines
- 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.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
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.
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 (3)
- [Section 4.1, Proposition 1 and Eq. (5)] 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 6.1.1 and Section 7] 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 5 and Section 6.2] 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.
minor comments (5)
- [Section 4.2 and Section 4.3] 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 6.1.1] 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 6.2] 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 8] 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 9 and Figure 3b] 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'.
Circularity Check
No significant circularity: the optimization objective directly targets the input moments, the reparametrization is proved self-contained, and self-citations are not load-bearing.
full rationale
The paper's central derivation is the differentiable reparametrization of PH distributions onto an unconstrained space (Section 4.1, Proposition 1). This is proved constructively: given any interior Markovian PH representation, the paper exhibits explicit pre-images a = log(alpha), gamma_i = sqrt(-T_ii), and Z = log(E), and verifies that the map pi recovers alpha and T. The same holds for the Coxian and Hyper-Erlang reparametrizations (Propositions 2 and 3). No step in this chain is defined in terms of the target moments, and no fitted parameter is renamed as a prediction. The moment-fitting objective (Eq. 6, 12, 18) is a weighted regression against the input moments themselves, so reporting success rates on those moments is a direct report of optimization quality, not a prediction of held-out quantities. The paper's evaluation uses moment sequences sampled from PH distributions (Section 6.1.1), which guarantees feasibility for sufficiently large n; this is an external-validity limitation for arbitrary real-world moment sequences, but it is not circular reasoning, because the optimizer receives only the moment targets and does not use the generating PH representation. Self-citations (e.g., Horvath and Telek 2024; Sherzer et al. 2024; A. Bobbio and Telek 2005) appear in the literature review and motivation, but none is invoked as the load-bearing justification for the reparametrization or the optimization results. The paper is self-contained against the benchmark of its own objective, and no circular step could be exhibited with a specific reduction of a result to its input.
Assumptions & free parameters
free parameters (4)
- Number of starting points s =
10000
- Early termination schedule =
2000 processes at step 500; 200 at step 5000; 20 at step 15000
- Moment weights w_i =
m_i^{-2}
- Shape fitting trade-off Q =
0.05
assumptions (4)
- standard math The i-th moment of a PH distribution equals i!(-1)^i alpha T^{-i} 1
- standard math The mapping pi is surjective onto the interior of the Markovian PH space
- domain assumption Boundary points of the PH space are not reachable but are not needed for numerical fitting
- domain assumption The non-singularity constraint on T can be omitted in numerical work
Cite this review
Pith. "Pith review of An Unconstrained Optimization Approach to Moment Fitting with Phase Type Distributions." pith.science (2026). https://pith.science/paper/O5IJNOGF
@misc{pith2026250520379,
author = {Pith},
title = {Pith review of: An Unconstrained Optimization Approach to Moment Fitting with Phase Type Distributions},
year = {2026},
howpublished = {\url{https://pith.science/paper/O5IJNOGF}},
note = {Machine review of arXiv:2505.20379}
}
read the original abstract
Phase type (PH) distributions are widely used in modeling and simulation due to their generality and analytical properties. In such settings, it is often necessary to construct a PH distribution that aligns with real-world data by matching a set of prescribed moments. Existing approaches provide either exact closed-form solutions or iterative procedures that may yield exact or approximate results. However, these methods are limited to matching a small number of moments using PH distributions with a small number of phases, or are restricted to narrow subclasses within the PH family. We address the problem of approximately fitting a larger set of given moments using potentially large PH distributions. We introduce an optimization methodology that relies on a re-parametrization of the Markovian representation, formulated in a space that enables unconstrained optimization of the moment-matching objective. This reformulation allows us to scale to significantly larger PH distributions and capture higher moments. Results on a large and diverse set of moment targets show that the proposed method is, in the vast majority of cases, capable of fitting as many as 20 moments to PH distributions with as many as 100 phases, with small relative errors on the order of under 0.5% from each target. We further demonstrate an application of the optimization framework where we search for a PH distribution that conforms not only to a given set of moments but also to a given shape. Finally, we illustrate the practical utility of this approach through a queueing application, presenting a case study that examines the influence of the i^{th} moment of the inter-arrival and service time distributions on the steady-state probabilities of the GI/GI/1 queue length.
Figures
Figures from the paper (7 more)
Reference graph
Works this paper leans on
-
[10]
On the canonical representation of phase type distributions.Perfor- mance Evaluation66, 8 (2009), 396–409.https://doi.org/10.1016/j.peva.2008.11.002Selected papers of the Fourth European Performance Engineering Workshop (EPEW) 2007 in Berlin. Mary A. Johnson and Michael R. Taaffe
-
[13]
A perishable inventory system with service facilities, MAP arrivals and PH — Service times.Journal of Systems Science and Systems Engineering16, 1 (01 Mar 2007), 62–73.https://doi.org/10.1007/s11518-006-5025-3 Raymond Marie
-
[17]
Approximating G(t)/GI/1 queues with deep learning.European Journal of Operational Research(2024).https://doi.org/10.1016/ j.ejor.2024.12.030 David Stanford and Wolfgang Fischer
work page 2024
-
[18]
The interdeparture-time distribution for each class in theP i Mi/Gi/1queue.Queueing Systems4, 3 (01 Jul 1989), 179–191.https://doi.org/10.1007/ BF02100265 M. Telek and G. Horváth
work page 1989
-
[19]
A minimal representation of Markov arrival processes and a moments matching method.Performance Evaluation64, 9 (2007), 1153–1168.https://doi.org/10.1016/j. peva.2007.06.001Performance
doi:10.1016/j 2007
-
[20]
In2005 International Conference on Dependable Systems and Networks (DSN’05)
A novel approach for fitting probability distributions to real trace data with the EM algorithm. In2005 International Conference on Dependable Systems and Networks (DSN’05). 712–721.https://doi.org/10.1109/DSN.2005.11 A. Thummler, P. Buchholz, and M. Telek
-
[21]
A Novel Approach for Phase-Type Fitting with the EM Algorithm.IEEE Transactions on Dependable and Secure Computing3, 3 (2006), 245–258.https: //doi.org/10.1109/TDSC.2006.27 Ward Whitt
-
[22]
Approximating a Point Process by a Renewal Process, I: Two Basic Meth- ods.Operations Research30, 1 (1982), 125–147.https://doi.org/10.1287/opre.30.1.125 arXiv:https://doi.org/10.1287/opre.30.1.125 19 Appendix A. Evaluation of previous methods In this section, we evaluate the optimization method proposed in Buchholz and Kriege (2009)3 as a function of a n...
Show all 22 references
-
[1975]
Approximate Analysis of Central Server Models.IBM Journal of Research and Development19, 3 (1975), 301–313.https://doi.org/10.1147/rd.193.0301 Eliran Sherzer, Opher Baron, Dmitry Krass, and Yehezkel Resheff
1975 doi
-
[1980]
Calculating equilibrium probabilities forλ(n)/Ck/1/N queues.SIGMETRICS Per- form. Eval. Rev.9, 2 (May 1980), 117–125.https://doi.org/10.1145/1009375.806155 Takayuki Osogami and Mor Harchol-Balter
1980
-
[1982]
On the canonical representation of homogeneous markov processes modelling failure - time distributions.Microelectronics Reliability22, 3 (1982), 583–602.https://doi.org/10.1016/ 0026-2714(82)90033-6 Peter I. Frazier
1982
-
[1985]
On the Phase-Type Approximations of General Distributions.IIE Transactions17, 2 (1985), 110–116.https://doi.org/10.1080/07408178508975280 arXiv:https://doi.org/10.1080/07408178508975280 Opher Baron, Dmitry Krass, Arik Senderovich, and Eliran Sherzer
1985 doi
-
[1989]
Stochastic Models5, 4 (1989), 711–743
Matching moments to phase distributions: Mixtures of erlang distributions of common order.Communications in Statistics. Stochastic Models5, 4 (1989), 711–743. https://doi.org/10.1080/15326348908807131arXiv:https://doi.org/10.1080/15326348908807131 Mary A. Johnson and Michael R...
1989 doi
-
[1991]
Sivakumar, and G
An investigation of phase-distribution moment-matching algorithms for use in queueing models.Queueing Systems8, 1 (01 Dec 1991), 129–147.https: //doi.org/10.1007/BF02412246 Paul Manuel, B. Sivakumar, and G. Arivarignan
1991 doi
-
[2005]
Matching Three Moments with Minimal Acyclic Phase Type Distributions.Stochastic Models21, 2-3 (2005), 303–326.https://doi.org/10.1081/STM-200056210 arXiv:https://doi.org/10.1081/STM-200056210 Tayfur Altiok
2005 doi
-
[2006]
Closed form solutions for mapping general distributions to quasi-minimal PH distributions.Performance Evaluation63, 6 (2006), 524–552.https://doi.org/ 10.1016/j.peva.2005.06.002ModellingTechniquesandToolsforComputerPerformanceEvaluation. C. H. Sauer and K. M. Chandy
2006 doi
-
[2007]
Matching More Than Three Moments with Acyclic Phase Type Dis- tributions.Stochastic Models23, 2 (2007), 167–194.https://doi.org/10.1080/15326340701300712 arXiv:https://doi.org/10.1080/15326340701300712 Gábor Horváth and Miklós Telek
2007 doi
-
[2008]
7461), Dario Bini, Beatrice Meini, Vaidyanathan Ramaswami, Marie- Ange Remiche, and Peter Taylor (Eds.)
Interarrival Times Characterization and Fitting forMarkovianTrafficAnalysis.InNumerical Methods for Structured Markov Chains (Dagstuhl Seminar Proceedings (DagSemProc), Vol. 7461), Dario Bini, Beatrice Meini, Vaidyanathan Ramaswami, Marie- Ange Remiche, and Peter Taylor (Eds.)...
-
[2009]
In2009 Sixth international conference on the quantitative evaluation of systems
A heuristic approach for fitting MAPs to moments and joint moments. In2009 Sixth international conference on the quantitative evaluation of systems. IEEE, 53–62. Peter Buchholz, Jan Kriege, and Iryna Felko. 2014.Input modeling with phase-type distributions and Markov models: t...
2014
-
[2012]
2024.Phase type distributions: theory and applications
Analysis of a continuous time SM[K]/PH[K]/1/FCFS queue: Age process, sojourn times, and queue lengths.Journal of Systems Science and Complexity25, 1 (01 Feb 2012), 133–155.https: //doi.org/10.1007/s11424-012-9138-0 András Horváth and Miklós Telek. 2024.Phase type distributions...
2012 doi
-
[2018]
arXiv:1807.02811 [stat.ML] Harchol-Balter, Mor
A Tutorial on Bayesian Optimization. arXiv:1807.02811 [stat.ML] Harchol-Balter, Mor. 2013.Performance Modeling and Design of Computer Systems: Queueing Theory in Action. Cambridge University Press. Qiming He
2013 arXiv
-
[2024]
1287/ijoc.2022.0263 17 Peter Buchholz and Jan Kriege
Supervised ML for Solving the GI/GI/1 Queue.INFORMS Journal on Computing36, 3 (2024), 766–786.https://doi.org/10. 1287/ijoc.2022.0263 17 Peter Buchholz and Jan Kriege
2024
Reviewed August 7, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.