Building Accurate Workload Models Using Markovian Arrival Processes
Building Accurate Workload Models Using Markovian Arrival Processes Giuliano Casale Department of Computing Imperial College London [email protected]
BCAM Seminar – November 15th, 2012
1/51
Outline
Main topics covered in this seminar: • Phase-type (PH) distributions • Moment matching • Markovian arrival processes (MAP)
2/51
1. PH DISTRIBUTIONS
3/51
Continuous-Time Markov Chain (CTMC) Notation • m states • λi,j ≥ 0: (exponential) transition rate from state i to j P • λi = m j=1 λi,j : total outgoing rate from state i • Infinitesimal generator matrix: −λ1 λ1,2 . . . λ1,m λ2,1 −λ2 . . . λ2,m Q= . .. .. , Q1 = 0, .. .. . . . λn,1
λn,2
1 1 1= .. .
. . . −λm
1
• π(t) ˙ = π(t)Q, initial value π(0) • π(t) = π(0)e Qt : state probability vector at time t • πi (t): probability of the CTMC being in state i at time t • π(0): initial state probability vector
4/51
Example 1.1: CTMC Transient Analysis Initial state: π(0) = [0.9, 0.0, 0.1, 0.0]
−5 4 1 0 4 −7 2 1 Q= 2 3 −5 0 2 0 0 −2
Transient analysis: π(t) = π(0)e Qt = π(0)
∞ X (Qt)k k=0
k!
Transient Probabilities
A Sample Path 1 0.8 state probabilities
current state
4
3
2
0.6 0.4 0.2
1 0
5
10
15
20
25 time
30
35
40
45
0 0
0.5
1
1.5 time
2
2.5
3
5/51
Phase-Type Distribution (PH) • Q: CTMC with m = n + 1 states, last is absorbing (λn+1 = 0)
−λ1 λ1,2 . . . λ1,n λ1,n+1
λ2,1 −λ2 . . . λ2,n λ2,n+1 T t .. .. . . .. .. Q= = . . . . 0 0 . λ n,1 λn,2 . . . −λn λn,n+1 0 0 ... 0 0
• T: PH subgenerator matrix, T1 < 0 • t = −T1: exit vector
• “No mass at zero” assumption: π(0) = [α, 0], α1 = 1 References: [Neu89]
6/51
Example 1.2: Absorbing State Probability Initial state: π(0) = [α, 0.0] = [0.9, 0.0, 0.1, 0.0] −5 4 1 0 −5 4 1 0 4 −7 2 1 Q= , T = 4 −7 2 , t = −T1 = 1 2 3 −5 0 2 3 −5 0 0 0 0 0 1 absorbing state probability
1
state probabilities
0.8 0.6 0.4 0.2 0 0
5
10 time
15
0.8 0.6 0.4 0.2 0 0
5
10
15
time
7/51
PH: Fundamental Concepts Basic idea: model F (t) = Pr(event occurs in X ≤ t time units) as the probability mass absorbed in t time units in Q. • Semantics: entering the absorbing state models the occurrence of an event, e.g., the arrival of a TCP packet, the completion of a job, a non-exponential state transition in a complex system • Understanding the absorption dynamics: ∞ ∞ P (Tt)k P Tk−1 t k t π(t) = π(0)e Qt = π(0) k=0 k! k=1 k! 0
1
• Using the definition t = −T1, we get π(t) = [αe Tt , πn+1 (t)],
πn+1 (t) = 1 − αe Tt 1
where πn+1 (t) is the probability mass absorbed in t time units.
8/51
PH: Fundamental Formulas
PH distribution: F (t) = Pr(event occurs in X ≤ t) = 1 − αe Tt 1 def
• PH representation: (α, T) • Probability Density function: f (t) = αe Tt (−T)1 R +∞ • (Power) Moments: E [X k ] = 0 t k f (t)dt = k!α(−T)−k 1 • (−T)−1 = [τi,j ] ≥ 0 • τi,j : mean time spent in state j if the PH starts in state i • Median/Percentiles: no simple form, determined numerically.
9/51
Example 1.3: a 3-state PH distribution
−(λ + µ) λ µ 0 0 −λ λ , t = 0 , λ ≥ 0. α = [1, 0, 0], T = 0 0 −λ λ • Case µ → +∞: E [X k ] = k!λ−k (Exponential, c 2 = 1). • Case µ = λ: E [X k ] = (k + 1)!λ−k (Erlang-2, c 2 = 1/2) (k+2)! −k 2 λ delivers c 2 > 1
• Case µ = 0: E [X k ] = • No choice of µ
(Erlang-3, c 2 = 1/3).
def
Remarks: Squared coefficient of variation: c 2 = Var [X ]/E [X ]2
10/51
PH: “Family Picture” - n ≤ 2 states c2 1
α [1]
1 2
[1, 0]
Hypo-Exp.
[ 12 , 1)
[1, 0]
Hyper-Exp.
[1, +∞)
[α1 , α2 ]
Coxian/APH
[ 12 , +∞)
[1, 0]
General
[ 12 , +∞)
[α1 , α2 ]
Exponential Erlang
T −λ −λ λ 0 −λ −λ1 λ1 −λ2 0 −λ1 0 −λ2 0 −λ1 p1 λ1 −λ2 0 −λ1 p1 λ1 p2 λ2 −λ2
Subset of Hyper-Exp. Hypo-Exp. Coxian/APH Coxian/APH General
def Remarks: α1 + α2 = 1, c 2 = Var [X ]/E [X ]2
11/51
PH: “Family Picture” - Examples — n = 3 states c2
α
Hyper-Erlang
[ 12 , +∞)
[α1 , α2 , 0]
Coxian/APH
[ 13 , +∞)
[1, 0, 0]
Circulant
[ 13 , +∞)
[α1 , α2 , α3 ]
General
[ 13 , +∞)
[α1 , α2 , α3 ]
T −λ1 0 0 −λ1 0 0 −λ1 0 p λ 31 3 −λ1 p21 λ2 p31 λ3
0 −λ2 0 p1 λ 1 −λ2 0 p12 λ1 −λ2 0 p12 λ1 −λ2 p32 λ3
0 λ2 −λ2 0 p2 λ2 −λ3 0 p23 λ2 −λ3 p13 λ1 p23 λ2 −λ3
def Remarks: α1 + α2 + α3 = 1, c 2 = Var [X ]/E [X ]2
12/51
Example 1.4: Reducing to Coxian Form Algorithms exist to reduce a PH to Coxian form. • With n = 2 states (PH(2) models) this can be done analytically −25 0 0 0 • Hyper-Exponential: α = [0.99, 0.01], T = , 0 −5 F 0 (t) = 1 − 0.99e −25t − 0.01e −5t −λ1 p1 λ1 • Coxian: α = [1, 0], T = 0 −λ2 • Symbolic analysis gives that in the 2-state Coxian F (t) = 1 − M1 e −λ1 t − (1 − M1 )e −λ2 t ,
def
M1 = 1 −
λ1 p1 λ1 − λ2
• Thus the two models are equivalent if λ1 = 25, λ2 = 5, and p1 = 0.008 such that M1 = 0.99.
13/51
Example 1.4: Reducing to Coxian Form Compare moments of Hyper-Exponential and Coxian: E [X ] E [X 2 ] E [X 3 ] E [X 4 ] .. .
Hyper-Exp 41.600 · 10−3 3.968 · 10−3 860.160 · 10−6 444.826 · 10−6 .. .
Coxian 41.600 · 10−3 3.968 · 10−3 860.160 · 10−6 444.826 · 10−6 .. .
Key message: differences between (α0 , T0 ) and (α, T) are deceptive! In general, (α, T) is a redundant representation. Redundancy problem: how many degrees of freedom in PH distributions? How to cope with redundant parameters?
14/51
PH: Degrees of Freedom • PH Moments: E [X k ] = k!α(−T)−k 1 • A of order n, characteristic polynomial: • φ(θ) = det(θI − A) • φ(A) = An + m1 An−1 + . . . + mn−1 A + mn I = 0 • A = (−T)−1 implies that PH moments are linearly-dependent E [X n−1 ] E [X n ] + m1 + . . . + mn−1 E [X ] + mn = 0 n! (n − 1)! • Thus, a PH offers up to 2n − 1 degrees of freedom (df) for fitting a workload distribution (n − 1 moments + n terms mj ). • n = 2 states ⇒ 3 df, n = 3 ⇒ 5 df, n = 4 ⇒ 7 df, ... References: [TelH07],[CasZS07]
15/51
PH: Algebra of Random Variables and Closure Properties • PH: the smallest family of distributions on 1 • Obtain m1 ,. . .,mn by solving the linear system e n] e n−1 ] E [X [X + m1 E(n−1)! + . . . + mn = 0 k! n+1 n e e E [X ] + m1 E [X ] + . . . + mn E [X e] = 0 (k+1)!
E [Xe 2n−1 ] (2k−1)!
.. .
n!
.. .
e 2n−2
e n−1
[X ] E [X ] + m1 E(2K −2)! + . . . + mn (n−1)! = 0
• Obtain θi as roots of φ(θ) = θn + m1 θn−1 + . . . + mn−1 θ + mn • Obtain Mi,1 from the spectral characterization given θi and e k ]. E [X −1 • Output: α = (M1,1 , . . . , MK ,1 ), T = −diag(θ1−1 , . . . , θK ) References: [CasZS08b]
23/51
Example 2.1: Prony’s Method Trace: LBL-TCP-3, http://ita.ee.lbl.gov MATLAB Code, PH(3): E=[1,5.3109e-003,89.7845e-006,3.0096e-006,163.1872e-009,12.9061e-009]; f=factorial(0:6); A=[E(4:-1:1)./f(4:-1:1);E(5:-1:2)./f(5:-1:2); E(6:-1:3)./f(6:-1:3);1,0,0,0]; b=[0;0;0;1]; m=(A\b)’, theta=roots(m)’ M=([theta*f(2);theta.ˆ2*f(3);theta.ˆ3*f(4)]\E(2:4)’)’
α = [M1,1 , M2,1 , M3,1 ] = [26.2074, 384.7137, 589.0789] · 10−3 ; T = −diag(θ1−1 , θ2−1 , θ3−1 ) = −diag(49.9628, 110.5089, 451.3770) BC−pOct89
LBL−TCP−3 0
−1
−1
−2 −3 −4 −5 −6 −5
Trace PH(3) −4
−3 −2 −1 Interarrival Time t [log]
−2 −3 −4 −5 −6
0
CCDF 1−F(t) [log]
0 −1 CCDF 1−F(t) [log]
CCDF 1−F(t) [log]
LBL−PKT−5 0
−7 −5
Trace PH(3) −4
−3 −2 −1 Interarrival Time t [log]
−2 −3 −4 −5 −6
0
−7 −5
Trace PH(3) −4
−3 −2 −1 Interarrival Time t [log]
0
24/51
PH: Approximate Moment Matching • Relatively few techniques for approximate moment matching • Kronecker product composition for PH (KPC-PH): • X1 ∼ (α, T) of order n, X2 ∼ (β, S) of order m def
def
• X = X1 ⊗ X2 ∼ (γ, R), γ = α ⊗ β, R = (−T) ⊗ S • (γ, R) is PH only if S = −diag(λ1 , . . . , λn ) • Divide-and-conquer approximate moment matching: E [X k ] =k!(α ⊗ β)(−((−T) ⊗ S))−k (1n ⊗ 1m ) =k!(α(−T)−k 1n )(β(−S)−k 1m ) =E [X1k ]E [X2k ]/k! References: [CasZS08]
25/51
Example 2.2: KPC-PH – Increased Degrees of Freedom • X1 : PH(2), E [X1 ] = 1, E [X12 ] = 10, E [X13 ] = 200 • X2 : PH(2), E [X2 ] = 1, E [X22 ] = 10, E [X23 ] = 3200 • Y1 : PH(2), E [Y1 ] = 1, E [Y12 ] = 3.2691, E [Y13 ] = 200 • Y2 : PH(2), E [Y2 ] = 1, E [Y22 ] = 30.589, E [Y23 ] = 3200 • Z1 : PH(2), E [Z1 ] = 1, E [Z12 ] = 20, E [Z13 ] = 200 • Z2 : PH(2), E [Z2 ] = 1, E [Z22 ] = 5, E [Z23 ] = 3200 E [X ] E [X 2 ] E [X 3 ] E [X 4 ] E [X 5 ] E [X 6 ] .. .
X1 ⊗ X2 1.0000 5.0000 · 101 1.0667 · 105 7.2362 · 108 5.2745 · 1012 4.8979 · 1016 .. .
Y1 ⊗ Y2 1.0000 5.0000 · 101 1.0667 · 105 7.2362 · 108 6.3130 · 1012 6.6129 · 1016 .. .
Z1 ⊗ Z2 1.0000 5.0000 · 101 1.0667 · 105 3.7813 · 108 1.6796 · 1012 8.9531 · 1015 .. .
26/51
PH: Generalized KPC-PH Technique Generalization: N N N • X = Jj=1 Xj ∼ ( Jj=1 αj , (−1)J−1 Jj=1 Tj ) • X is PH if J − 1 subgenerators are diagonal Q E [X k ] • E [X k ] = k! Jj=1 k!j KPC-Toolbox:
http://www.cs.wm.edu/MAPQN/kpctoolbox.html
• Support for exact and approximate moment matching • Determine (αj , Tj ) by numerical optimization • Search on moments directly instead of (α, T) parameters References: [CasZS08]
27/51
Example 2.3: KPC-Toolbox Fitting • Internet Traffic Archive Trace: BC-Aug-p89,
http://ita.ee.lbl.gov
• Prony’s method fails if n > 3: θ2 = 0.0026 + i0.0179 • How to find a perturbation of the moment sets that delivers more accurate result? 0
10
KPC-PH, PH(4): 39.630, 155.67) α =[0.66082, 0.31138,
CCDF 1−F(t) [log]
T = − diag(283.677, 1114.3,
−2
10
−4
10
Trace PH(3) − Prony PH(4) − KPC−PH
−6
0.01889, 0.00890]
10
−8
10 −5 10
−4
10
−3
−2
10 10 Interarrival Time [log]
−1
10
0
10
28/51
PH: Other Publicly Available Tools for PH Fitting EMpht (1996) —
http://home.imf.au.dk/asmus/pspapers.html:
• EM algorithm for ML fitting, based on Runge-Kutta methods • Local optimization technique jPhase (2006) —
http://copa.uniandes.edu.co/software/jmarkov/index.html:
• Java library — ML and canonical form fitting algorithms PhFit (2002) —
http://webspn.hit.bme.hu/∼telek/tools.htm:
• Separate fit of distribution body and tail • Both continuous and discrete ML distributions G-FIT (2007) —
http://ls4-www.cs.uni-dortmund.de/home/thummler/gfit.tgz:
• Hyper-Erlang PHs used as building block • Automatic aggregation of large traces, dramatic speed-up of computational times compared to EMpht
29/51
3. MARKOVIAN ARRIVAL PROCESS
30/51
Time Series Analysis
Notation: def
• {t0 = 0, t1 , t2 , t3 , . . .} : sequence of arrival times of events in the real system def
• Xk = tk − tk−1 : inter-arrival time between arrival of the (k − 1)-th and the k-th events. • tk and Xk may not be directly observable, e.g., aggregate data
31/51
Stochastic Process Descriptions
• Inter-arrival process: models sequence of values Xk • Natural description for unaggregated traces • Enables reasoning on individual events (e.g., response time distributions, covariance of successive arrivals, ...) • Counting process: models number of arrivals in interval [0, t] • Preferred for aggregate data (e.g., packet counts) • Enables reasoning on the volumes of events at timescale t • The two descriptions are equivalent in theory; in practice they carry independent information when fitted on a dataset
32/51
Sequence of PH Inter-Arrival Times
• PH-Renewal process: • for k = 1, 2, 3, . . . • initialize a PH with subgenerator T according to α • generate Xk as the time to absorption in (α, T) • tk = tk−1 + Xk is the arrival time of the k-th event • Limitation: every time the PH is restarted independently of the past. No way to define time-varying patterns, e.g., periodicities, burstiness, ... References: [Neu89]
33/51
PH-R: Counting Process
• PH-Renewal Process State (N(t), J(t)) • N(t): event counter increased upon arrival/restart events • J(t): PH state within current level α T t·α 0 0 0 0 T t · α 0 π c (0) = 0 Qc = 0 0 T t·α .. .. .. .. .. . . . . .
. . . → N(t) = 0 . . . → N(t) = 1 . . . → N(t) = 2 .. .. . .
34/51
Markovian Arrival Process (MAP) • MAP: generalizes the PH-Renewal construction by considering restarts that depend on the exit state of the previous sample. • A technical device to introduce “memory” in the time series. • MAP Counting Process: D0 D1 0 0 . . . → N(t) = 0 α 0 D0 D1 0 . . . → N(t) = 1 0 π c (0) = 0 Qc = 0 0 D0 D1 . . . → N(t) = 2 .. .. .. .. . . . . .. . . . . . . . • D0 : same as T, transitions do not increase N(t) • D1 : generalizes t · α, transitions increase N(t) by 1 • Batch MAP (BMAP): Db , b > 1, increasing N(t) by b units References: [Neu89]
35/51
Example 3.1: Sequence of MAP(2) Inter-Arrival Times −10 3 5 D0 = , D1 = 5 −5 0
2 0
• k = 1; start in a random state according to α, e.g., state 2. • Xk = 0 • Xk = Xk + r , r ∼ exp(5). Jump to state 1. • Xk = Xk + r , r ∼ exp(10). • u ∼ uniform(0, 1). 3 • if u ∈ [0, 10 ) jump to state 2 3 3+5 • if u ∈ [ 10 , 10 ): save Xk ; Xk+1 = 0; restart from state 1. 3+5+2 • if u ∈ [ 3+5 10 , 10 ]: save Xk ; Xk+1 = 0; restart from state 2
• ...
36/51
Example 3.2: Temporal Dependence in MAPs
−1 D0 = 0 0
0 −10 0
0 p 1−p 0 , D1 = 0 10p −100 100(1 − p) 0
0 10(1 − p) 100p
p = 0.90
p = 0.99
p = 0.10
2
2
10
10
0
−2
10
−4
10
−2
10
−4
10
−6
0
0
10
sample value − Sk
sample value − Sk
sample value − Sk
0
10
10
2
10
400 600 sample number − k
800
1000
10
0
−2
10
−4
10
−6
200
10
−6
200
400 600 sample number − k
800
1000
10
0
200
400 600 sample number − k
800
1000
37/51
MAP: “Family Picture” - n ≤ 2 states Name Poisson Erlang Renewal Process Hyper-exp. Renewal Process Interrupted Poisson Process MMPP Acyclic MAP(2) MAP(2)
D0 −λ −λ λ 0 −λ −λ1 0 0 −λ2 −λ1 λ1,2 λ −λ2 2,1 −λ1 λ1,2 λ −λ2 2,1 −λ1 λ1,2 0 −λ2 −λ1 λ1,2 λ2,1 −λ2
D1 λ 0 0 λ 0 pλ1 qλ1 pλ2 qλ2 ∗ λ1,1 λ∗1,2 0 0 ∗ λ1,1 0 0 λ∗2,2 ∗ λ1,1 0 λ∗ λ∗2,2 ∗2,1 λ1,1 λ∗1,2 λ∗2,1 λ∗2,2
Remarks: p + q = 1
38/51
MAP: “Family Picture” - Examples – n = 3 states Name −λ1 0 Hyper-Exp. Renewal 0 −λ1 Circulant MMPP(3) λ2,1 0 −λ1 λ2,1 MMPP(3) λ3,1 −λ1 0 Hyper-Exp. MAP 0 −λ1 λ2,1 MAP(3) λ3,1 Remarks: p + q + r = 1
D0 0 −λ2 0 0 −λ2 λ3,2 λ1,2 −λ2 λ3,2 0 −λ2 0 λ1,2 −λ2 λ3,2
0 0 −λ3 λ1,1 0 −λ3 λ1,3 λ2,3 −λ3 0 0 −λ3 λ1,3 λ2,3 −λ3
λ1 p λ1 p λ p ∗3 λ1,1 0 0 ∗ λ1,1 0 0 ∗ λ1,1 λ∗2,1 λ∗ 3,1 λ∗1,1 λ∗2,1 λ∗3,1
D1 λ1 q λ2 q λ3 q 0 λ∗2,2 0 0 λ∗2,2 0 λ∗1,2 λ∗2,2 λ∗3,2 λ∗1,2 λ∗2,2 λ∗3,2
λ1 r λ2 r λ3 r 0 0 λ∗3,3 0 0 λ∗3,3 λ∗1,3 λ∗2,3 λ∗3,3 λ∗1,3 λ∗2,3 λ∗3,3
39/51
MAP: Key Formulas for Inter-Arrival Times
• MAP representation: (D0 , D1 ) • Embedded chain: P = (−D0 )−1 D1 • Interval-stationary initial vector: α = αP Distribution of samples: • (Marginal) Distribution: F (t) = 1 − αe D0 t 1
• (Marginal) Density: f (t) = αe D0 t (−D0 )1 = αe D0 t D1 1
• (Marginal) Moments: E [X k ] = k!α(−D0 )−k 1 • Degrees of freedom of distribution: 2n − 1
40/51
MAP: Key Formulas for Inter-Arrival Times Sequence of samples: • Joint Density Function: Pr(X1 = x1 , X2 = x2 , . . . , Xq = xq ) = αe D0 x1 D1 e D0 x2 D1 · · · e D0 xq D1 1
• D1 → D1 Ph−1 for samples spaced by h lags • Joint Moments: E [X1k1 X2k2 · · · Xq q ] = k1 ! · · · kq !α(−D0 )−k1 P(−D0 )−k2 · · · (−D0 )−kq 1 k
• P → Ph for samples spaced by h lags • Analysis often limited to second-order properties. • Autocorrelation function: ρh =
α(−D0 )−1 Ph (−D0 )−1 1 − α(−D0 )−1 1 E [X1 X1+h ] − E [X ]2 = E [X 2 ] − E [X ]2 2α(−D0 )−2 1 − α(−D0 )−1 1 References: [TelH07]
41/51
4. INTER-ARRIVAL PROCESS FITTING
Includes joint work with Evgenia Smirni (IBM Research, William & Mary)
42/51
MAP: Approximate Inter-Arrival Process Fitting KPC readily generalizes to MAP random variables N (j) N (j) N • X = Jj Xj ∼ ((−1)J−1 Jj D0 , Jj D1 ) (j)
• X is a MAP if J − 1 D0 are diagonal (hyper-exp. moments) Moments and Joint Moments: Q E [X k ] • E [X k ] = k! Jj=1 k!j u QJ E [Xj,ik Xj,i+h ] u ] = k!u! • E [Xik Xi+h , ∀ lags i, h j=1 k!u! u v QJ E [Xj,ik Xj,i+h Xj,i+h+l ] u Xv • E [Xik Xi+h , ∀ lags i, h, l j=1 i+h+l ] = k!u!v ! k!u!v ! • ... References: [CasZS07]
43/51
MAP: KPC-Toolbox Autocorrelation decay rates and embedded chain: Q N • γi = Jj=1 γj,i , P = Jj=1 Pj Second-order properties follow recursively from those of X ⊗ Y : 1 + c 2 = 21 (1 + cX2 )(1 + cY2 ) 2 2 2 2 cX cY cX cY X Y Y ρk = ρk + ρk + ρX k ρk 2 2 c c c2 KPC-Toolbox – MAP Fitting: • Optimization-based second-order and third-order fitting • Xj ∼ MAP(2), known feasibility region for c 2 and ρk • Fitting of c 2 and ρk disjoint from first-order and third-order • Residual df spent to fit third-order moments E [Xi Xi+h Xi+h+l ] References: [CasZS08]
44/51
Example 4.1: Second-Order vs Third-Order Fitting • Second-Order: approximately match E [Xi Xi+h ] (equiv. ρh ) • Third-Order: approx. match E [Xi Xi+h ] and E [Xi Xi+h Xi+h+l ] LBL−PKT−5
LBL−PKT−5
Pr(queue ≥ x) [log]
Trace Autocorrelation Identical Fittings
−1
−2
k
ρ IAT autocorrelations [log]
0
Trace Simulation Second−Order Fit Third−Order Fit
0
−2
−4
−3 0
1
2 3 lag k [log]
4
0
2 x [log]
4
References: [AndN02],[CasZS07]
45/51
Example 4.2: KPC MAP Fitting Comparison of inter-arrival and counting process methods: • KPC: second-order and third-order fit of inter-arrival process • Superposition: fit Hurst coefficient for counts SWeb − MAP/D/1 − Utilization 80%
SWeb − Autocorrelation Fitting
0
0
−1 Pr(queue ≥ x) [log]
Autocorrelation ρk [log]
Web MAP(16) kpc −1
−2
−2 −3 −4 −5
−3 0
1
2 3 Lag k [log]
4
5
−6 0
Trace SUP KPC 2
4
6
x [log]
References: [AndN98],[CasZS07]
46/51
CONCLUSION
47/51
Summary
• PHs and MAPs are tractable models for characterizing empirical data using Markov chains • PHs and MAPs are compositional with other Markov models, e.g., queueing systems • PHs and MAPs closed under several operations (mixture, convolution, kronecker sum, kronecker product, ...) • Models with n ≤ 3 states are analytically tractable • Larger models can be dealt with using divide-and-conquer methods (KPC)
48/51
REFERENCES
49/51
References and Additional Readings 1/2 [AldS87] D. Aldous and L. Shepp. The least variable phase type distribution is Erlang. Comm. Stat.-Stochastic Models, 3(3):467–473, 1987. [AndN98] A. T. Andersen and B. F. Nielsen. A Markovian approach for modeling packet traffic with long-range dependence. IEEE JSAC, 16(5):719–732, 1998. [AndN02] A. T. Andersen and B. F. Nielsen. On the use of second-order descriptors to predict queueing behavior of MAPs. Naval Res. Log., 49(4):391–409, 2002. [BauBK10] Falko Bause, Peter Buchholz, Jan Kriege. ProFiDo - The Processes Fitting Toolkit Dortmund. QEST 2010: 87-96. [BuchK09] P. Buchholz, Jan Kriege. A Heuristic Approach for Fitting MAPs to Moments and Joint Moments. In Proc of QEST, 2009: 53-62, 2009. [BodHGT08] L. Bodrog and A. Heindl and G. Horv´ ath and M. Telek. A Markovian Canonical Form of Second-Order Matrix-Exponential Processes. European Journal of Operations Research, 190(457-477), 2008. [Bre78] J. Brewer. Kronecker products and matrix calculus in system theory. IEEE Trans. on Circuits and Systems, 25(9):772–781, 1978. [CasZS07] G. Casale, E. Zhang, and E. Smirni. Characterization of moments and autocorrelation in MAPs. ACM Perf. Eval. Rev., 35(1):27–29, 2007. [CasZS08] G. Casale, E. Zhang, and E. Smirni. KPC-toolbox: Simple yet effective trace fitting using markovian arrival processes. In Proc. of QEST, 83–92, 2008. [CasZS08b] G. Casale, E. Z. Zhang, and E. Smirni. Interarrival times characterization and fitting for Markovian traffic analysis. Technical Report WM-CS-2008-02, College of William and Mary, 2008.
50/51
References and Additional Readings 2/2 [HeiHG06] A. Heindl, G. Horvath, and K. Gross. Explicit Inverse Characterizations of Acyclic MAPs of Second Order. EPEW workshop, 108-122, 2006. [HeiML06] A. Heindl, K. Mitchell, and A. van de Liefvoort. Correlation bounds for second-order MAPs with application to queueing network decomposition. Performance Evaluation, 63(6):553–577, 2006. [HorT09] G. Horvath, M. Telek. On the canonical representation of phase type distributions. Performance Evaluation, 66(8):396–409, (2009). [MaiO92] R. S. Maier and C. A O’Cinneide. A closure characterisation of phase-type distributions. J. App. Probab., 29:92–103, 1992. [Neu89] M. F. Neuts. Structured Stochastic Matrices of M/G/1 Type and Their Applications. Marcel Dekker, New York, 1989. [Nie98] B. F. Nielsen. Note on the Markovian Arrival Process, Nov 1998. http://www2.imm.dtu.dk/courses/04441/map.pdf [OsoH06] T. Osogami, M. Harchol-Balter. Closed form solutions for mapping general distributions to quasi-minimal PH distributions. Performance Evaluation, 63(6): 524-552, 2006. [TelH03] M. Telek and A. Heindl. Matching moments for acyclic discrete and continuous phase-type distributions of second order. International Journal of Simulation, 3:47–57, 2003. [TelH07] M. Telek and G. Horv´ ath. A minimal representation of Markov arrival processes and a moments matching method. Performance Evaluation, 64(9-12):1153-1168, 2007.
51/51
BCAM Seminar – November 15th, 2012
1/51
Outline
Main topics covered in this seminar: • Phase-type (PH) distributions • Moment matching • Markovian arrival processes (MAP)
2/51
1. PH DISTRIBUTIONS
3/51
Continuous-Time Markov Chain (CTMC) Notation • m states • λi,j ≥ 0: (exponential) transition rate from state i to j P • λi = m j=1 λi,j : total outgoing rate from state i • Infinitesimal generator matrix: −λ1 λ1,2 . . . λ1,m λ2,1 −λ2 . . . λ2,m Q= . .. .. , Q1 = 0, .. .. . . . λn,1
λn,2
1 1 1= .. .
. . . −λm
1
• π(t) ˙ = π(t)Q, initial value π(0) • π(t) = π(0)e Qt : state probability vector at time t • πi (t): probability of the CTMC being in state i at time t • π(0): initial state probability vector
4/51
Example 1.1: CTMC Transient Analysis Initial state: π(0) = [0.9, 0.0, 0.1, 0.0]
−5 4 1 0 4 −7 2 1 Q= 2 3 −5 0 2 0 0 −2
Transient analysis: π(t) = π(0)e Qt = π(0)
∞ X (Qt)k k=0
k!
Transient Probabilities
A Sample Path 1 0.8 state probabilities
current state
4
3
2
0.6 0.4 0.2
1 0
5
10
15
20
25 time
30
35
40
45
0 0
0.5
1
1.5 time
2
2.5
3
5/51
Phase-Type Distribution (PH) • Q: CTMC with m = n + 1 states, last is absorbing (λn+1 = 0)
−λ1 λ1,2 . . . λ1,n λ1,n+1
λ2,1 −λ2 . . . λ2,n λ2,n+1 T t .. .. . . .. .. Q= = . . . . 0 0 . λ n,1 λn,2 . . . −λn λn,n+1 0 0 ... 0 0
• T: PH subgenerator matrix, T1 < 0 • t = −T1: exit vector
• “No mass at zero” assumption: π(0) = [α, 0], α1 = 1 References: [Neu89]
6/51
Example 1.2: Absorbing State Probability Initial state: π(0) = [α, 0.0] = [0.9, 0.0, 0.1, 0.0] −5 4 1 0 −5 4 1 0 4 −7 2 1 Q= , T = 4 −7 2 , t = −T1 = 1 2 3 −5 0 2 3 −5 0 0 0 0 0 1 absorbing state probability
1
state probabilities
0.8 0.6 0.4 0.2 0 0
5
10 time
15
0.8 0.6 0.4 0.2 0 0
5
10
15
time
7/51
PH: Fundamental Concepts Basic idea: model F (t) = Pr(event occurs in X ≤ t time units) as the probability mass absorbed in t time units in Q. • Semantics: entering the absorbing state models the occurrence of an event, e.g., the arrival of a TCP packet, the completion of a job, a non-exponential state transition in a complex system • Understanding the absorption dynamics: ∞ ∞ P (Tt)k P Tk−1 t k t π(t) = π(0)e Qt = π(0) k=0 k! k=1 k! 0
1
• Using the definition t = −T1, we get π(t) = [αe Tt , πn+1 (t)],
πn+1 (t) = 1 − αe Tt 1
where πn+1 (t) is the probability mass absorbed in t time units.
8/51
PH: Fundamental Formulas
PH distribution: F (t) = Pr(event occurs in X ≤ t) = 1 − αe Tt 1 def
• PH representation: (α, T) • Probability Density function: f (t) = αe Tt (−T)1 R +∞ • (Power) Moments: E [X k ] = 0 t k f (t)dt = k!α(−T)−k 1 • (−T)−1 = [τi,j ] ≥ 0 • τi,j : mean time spent in state j if the PH starts in state i • Median/Percentiles: no simple form, determined numerically.
9/51
Example 1.3: a 3-state PH distribution
−(λ + µ) λ µ 0 0 −λ λ , t = 0 , λ ≥ 0. α = [1, 0, 0], T = 0 0 −λ λ • Case µ → +∞: E [X k ] = k!λ−k (Exponential, c 2 = 1). • Case µ = λ: E [X k ] = (k + 1)!λ−k (Erlang-2, c 2 = 1/2) (k+2)! −k 2 λ delivers c 2 > 1
• Case µ = 0: E [X k ] = • No choice of µ
(Erlang-3, c 2 = 1/3).
def
Remarks: Squared coefficient of variation: c 2 = Var [X ]/E [X ]2
10/51
PH: “Family Picture” - n ≤ 2 states c2 1
α [1]
1 2
[1, 0]
Hypo-Exp.
[ 12 , 1)
[1, 0]
Hyper-Exp.
[1, +∞)
[α1 , α2 ]
Coxian/APH
[ 12 , +∞)
[1, 0]
General
[ 12 , +∞)
[α1 , α2 ]
Exponential Erlang
T −λ −λ λ 0 −λ −λ1 λ1 −λ2 0 −λ1 0 −λ2 0 −λ1 p1 λ1 −λ2 0 −λ1 p1 λ1 p2 λ2 −λ2
Subset of Hyper-Exp. Hypo-Exp. Coxian/APH Coxian/APH General
def Remarks: α1 + α2 = 1, c 2 = Var [X ]/E [X ]2
11/51
PH: “Family Picture” - Examples — n = 3 states c2
α
Hyper-Erlang
[ 12 , +∞)
[α1 , α2 , 0]
Coxian/APH
[ 13 , +∞)
[1, 0, 0]
Circulant
[ 13 , +∞)
[α1 , α2 , α3 ]
General
[ 13 , +∞)
[α1 , α2 , α3 ]
T −λ1 0 0 −λ1 0 0 −λ1 0 p λ 31 3 −λ1 p21 λ2 p31 λ3
0 −λ2 0 p1 λ 1 −λ2 0 p12 λ1 −λ2 0 p12 λ1 −λ2 p32 λ3
0 λ2 −λ2 0 p2 λ2 −λ3 0 p23 λ2 −λ3 p13 λ1 p23 λ2 −λ3
def Remarks: α1 + α2 + α3 = 1, c 2 = Var [X ]/E [X ]2
12/51
Example 1.4: Reducing to Coxian Form Algorithms exist to reduce a PH to Coxian form. • With n = 2 states (PH(2) models) this can be done analytically −25 0 0 0 • Hyper-Exponential: α = [0.99, 0.01], T = , 0 −5 F 0 (t) = 1 − 0.99e −25t − 0.01e −5t −λ1 p1 λ1 • Coxian: α = [1, 0], T = 0 −λ2 • Symbolic analysis gives that in the 2-state Coxian F (t) = 1 − M1 e −λ1 t − (1 − M1 )e −λ2 t ,
def
M1 = 1 −
λ1 p1 λ1 − λ2
• Thus the two models are equivalent if λ1 = 25, λ2 = 5, and p1 = 0.008 such that M1 = 0.99.
13/51
Example 1.4: Reducing to Coxian Form Compare moments of Hyper-Exponential and Coxian: E [X ] E [X 2 ] E [X 3 ] E [X 4 ] .. .
Hyper-Exp 41.600 · 10−3 3.968 · 10−3 860.160 · 10−6 444.826 · 10−6 .. .
Coxian 41.600 · 10−3 3.968 · 10−3 860.160 · 10−6 444.826 · 10−6 .. .
Key message: differences between (α0 , T0 ) and (α, T) are deceptive! In general, (α, T) is a redundant representation. Redundancy problem: how many degrees of freedom in PH distributions? How to cope with redundant parameters?
14/51
PH: Degrees of Freedom • PH Moments: E [X k ] = k!α(−T)−k 1 • A of order n, characteristic polynomial: • φ(θ) = det(θI − A) • φ(A) = An + m1 An−1 + . . . + mn−1 A + mn I = 0 • A = (−T)−1 implies that PH moments are linearly-dependent E [X n−1 ] E [X n ] + m1 + . . . + mn−1 E [X ] + mn = 0 n! (n − 1)! • Thus, a PH offers up to 2n − 1 degrees of freedom (df) for fitting a workload distribution (n − 1 moments + n terms mj ). • n = 2 states ⇒ 3 df, n = 3 ⇒ 5 df, n = 4 ⇒ 7 df, ... References: [TelH07],[CasZS07]
15/51
PH: Algebra of Random Variables and Closure Properties • PH: the smallest family of distributions on 1 • Obtain m1 ,. . .,mn by solving the linear system e n] e n−1 ] E [X [X + m1 E(n−1)! + . . . + mn = 0 k! n+1 n e e E [X ] + m1 E [X ] + . . . + mn E [X e] = 0 (k+1)!
E [Xe 2n−1 ] (2k−1)!
.. .
n!
.. .
e 2n−2
e n−1
[X ] E [X ] + m1 E(2K −2)! + . . . + mn (n−1)! = 0
• Obtain θi as roots of φ(θ) = θn + m1 θn−1 + . . . + mn−1 θ + mn • Obtain Mi,1 from the spectral characterization given θi and e k ]. E [X −1 • Output: α = (M1,1 , . . . , MK ,1 ), T = −diag(θ1−1 , . . . , θK ) References: [CasZS08b]
23/51
Example 2.1: Prony’s Method Trace: LBL-TCP-3, http://ita.ee.lbl.gov MATLAB Code, PH(3): E=[1,5.3109e-003,89.7845e-006,3.0096e-006,163.1872e-009,12.9061e-009]; f=factorial(0:6); A=[E(4:-1:1)./f(4:-1:1);E(5:-1:2)./f(5:-1:2); E(6:-1:3)./f(6:-1:3);1,0,0,0]; b=[0;0;0;1]; m=(A\b)’, theta=roots(m)’ M=([theta*f(2);theta.ˆ2*f(3);theta.ˆ3*f(4)]\E(2:4)’)’
α = [M1,1 , M2,1 , M3,1 ] = [26.2074, 384.7137, 589.0789] · 10−3 ; T = −diag(θ1−1 , θ2−1 , θ3−1 ) = −diag(49.9628, 110.5089, 451.3770) BC−pOct89
LBL−TCP−3 0
−1
−1
−2 −3 −4 −5 −6 −5
Trace PH(3) −4
−3 −2 −1 Interarrival Time t [log]
−2 −3 −4 −5 −6
0
CCDF 1−F(t) [log]
0 −1 CCDF 1−F(t) [log]
CCDF 1−F(t) [log]
LBL−PKT−5 0
−7 −5
Trace PH(3) −4
−3 −2 −1 Interarrival Time t [log]
−2 −3 −4 −5 −6
0
−7 −5
Trace PH(3) −4
−3 −2 −1 Interarrival Time t [log]
0
24/51
PH: Approximate Moment Matching • Relatively few techniques for approximate moment matching • Kronecker product composition for PH (KPC-PH): • X1 ∼ (α, T) of order n, X2 ∼ (β, S) of order m def
def
• X = X1 ⊗ X2 ∼ (γ, R), γ = α ⊗ β, R = (−T) ⊗ S • (γ, R) is PH only if S = −diag(λ1 , . . . , λn ) • Divide-and-conquer approximate moment matching: E [X k ] =k!(α ⊗ β)(−((−T) ⊗ S))−k (1n ⊗ 1m ) =k!(α(−T)−k 1n )(β(−S)−k 1m ) =E [X1k ]E [X2k ]/k! References: [CasZS08]
25/51
Example 2.2: KPC-PH – Increased Degrees of Freedom • X1 : PH(2), E [X1 ] = 1, E [X12 ] = 10, E [X13 ] = 200 • X2 : PH(2), E [X2 ] = 1, E [X22 ] = 10, E [X23 ] = 3200 • Y1 : PH(2), E [Y1 ] = 1, E [Y12 ] = 3.2691, E [Y13 ] = 200 • Y2 : PH(2), E [Y2 ] = 1, E [Y22 ] = 30.589, E [Y23 ] = 3200 • Z1 : PH(2), E [Z1 ] = 1, E [Z12 ] = 20, E [Z13 ] = 200 • Z2 : PH(2), E [Z2 ] = 1, E [Z22 ] = 5, E [Z23 ] = 3200 E [X ] E [X 2 ] E [X 3 ] E [X 4 ] E [X 5 ] E [X 6 ] .. .
X1 ⊗ X2 1.0000 5.0000 · 101 1.0667 · 105 7.2362 · 108 5.2745 · 1012 4.8979 · 1016 .. .
Y1 ⊗ Y2 1.0000 5.0000 · 101 1.0667 · 105 7.2362 · 108 6.3130 · 1012 6.6129 · 1016 .. .
Z1 ⊗ Z2 1.0000 5.0000 · 101 1.0667 · 105 3.7813 · 108 1.6796 · 1012 8.9531 · 1015 .. .
26/51
PH: Generalized KPC-PH Technique Generalization: N N N • X = Jj=1 Xj ∼ ( Jj=1 αj , (−1)J−1 Jj=1 Tj ) • X is PH if J − 1 subgenerators are diagonal Q E [X k ] • E [X k ] = k! Jj=1 k!j KPC-Toolbox:
http://www.cs.wm.edu/MAPQN/kpctoolbox.html
• Support for exact and approximate moment matching • Determine (αj , Tj ) by numerical optimization • Search on moments directly instead of (α, T) parameters References: [CasZS08]
27/51
Example 2.3: KPC-Toolbox Fitting • Internet Traffic Archive Trace: BC-Aug-p89,
http://ita.ee.lbl.gov
• Prony’s method fails if n > 3: θ2 = 0.0026 + i0.0179 • How to find a perturbation of the moment sets that delivers more accurate result? 0
10
KPC-PH, PH(4): 39.630, 155.67) α =[0.66082, 0.31138,
CCDF 1−F(t) [log]
T = − diag(283.677, 1114.3,
−2
10
−4
10
Trace PH(3) − Prony PH(4) − KPC−PH
−6
0.01889, 0.00890]
10
−8
10 −5 10
−4
10
−3
−2
10 10 Interarrival Time [log]
−1
10
0
10
28/51
PH: Other Publicly Available Tools for PH Fitting EMpht (1996) —
http://home.imf.au.dk/asmus/pspapers.html:
• EM algorithm for ML fitting, based on Runge-Kutta methods • Local optimization technique jPhase (2006) —
http://copa.uniandes.edu.co/software/jmarkov/index.html:
• Java library — ML and canonical form fitting algorithms PhFit (2002) —
http://webspn.hit.bme.hu/∼telek/tools.htm:
• Separate fit of distribution body and tail • Both continuous and discrete ML distributions G-FIT (2007) —
http://ls4-www.cs.uni-dortmund.de/home/thummler/gfit.tgz:
• Hyper-Erlang PHs used as building block • Automatic aggregation of large traces, dramatic speed-up of computational times compared to EMpht
29/51
3. MARKOVIAN ARRIVAL PROCESS
30/51
Time Series Analysis
Notation: def
• {t0 = 0, t1 , t2 , t3 , . . .} : sequence of arrival times of events in the real system def
• Xk = tk − tk−1 : inter-arrival time between arrival of the (k − 1)-th and the k-th events. • tk and Xk may not be directly observable, e.g., aggregate data
31/51
Stochastic Process Descriptions
• Inter-arrival process: models sequence of values Xk • Natural description for unaggregated traces • Enables reasoning on individual events (e.g., response time distributions, covariance of successive arrivals, ...) • Counting process: models number of arrivals in interval [0, t] • Preferred for aggregate data (e.g., packet counts) • Enables reasoning on the volumes of events at timescale t • The two descriptions are equivalent in theory; in practice they carry independent information when fitted on a dataset
32/51
Sequence of PH Inter-Arrival Times
• PH-Renewal process: • for k = 1, 2, 3, . . . • initialize a PH with subgenerator T according to α • generate Xk as the time to absorption in (α, T) • tk = tk−1 + Xk is the arrival time of the k-th event • Limitation: every time the PH is restarted independently of the past. No way to define time-varying patterns, e.g., periodicities, burstiness, ... References: [Neu89]
33/51
PH-R: Counting Process
• PH-Renewal Process State (N(t), J(t)) • N(t): event counter increased upon arrival/restart events • J(t): PH state within current level α T t·α 0 0 0 0 T t · α 0 π c (0) = 0 Qc = 0 0 T t·α .. .. .. .. .. . . . . .
. . . → N(t) = 0 . . . → N(t) = 1 . . . → N(t) = 2 .. .. . .
34/51
Markovian Arrival Process (MAP) • MAP: generalizes the PH-Renewal construction by considering restarts that depend on the exit state of the previous sample. • A technical device to introduce “memory” in the time series. • MAP Counting Process: D0 D1 0 0 . . . → N(t) = 0 α 0 D0 D1 0 . . . → N(t) = 1 0 π c (0) = 0 Qc = 0 0 D0 D1 . . . → N(t) = 2 .. .. .. .. . . . . .. . . . . . . . • D0 : same as T, transitions do not increase N(t) • D1 : generalizes t · α, transitions increase N(t) by 1 • Batch MAP (BMAP): Db , b > 1, increasing N(t) by b units References: [Neu89]
35/51
Example 3.1: Sequence of MAP(2) Inter-Arrival Times −10 3 5 D0 = , D1 = 5 −5 0
2 0
• k = 1; start in a random state according to α, e.g., state 2. • Xk = 0 • Xk = Xk + r , r ∼ exp(5). Jump to state 1. • Xk = Xk + r , r ∼ exp(10). • u ∼ uniform(0, 1). 3 • if u ∈ [0, 10 ) jump to state 2 3 3+5 • if u ∈ [ 10 , 10 ): save Xk ; Xk+1 = 0; restart from state 1. 3+5+2 • if u ∈ [ 3+5 10 , 10 ]: save Xk ; Xk+1 = 0; restart from state 2
• ...
36/51
Example 3.2: Temporal Dependence in MAPs
−1 D0 = 0 0
0 −10 0
0 p 1−p 0 , D1 = 0 10p −100 100(1 − p) 0
0 10(1 − p) 100p
p = 0.90
p = 0.99
p = 0.10
2
2
10
10
0
−2
10
−4
10
−2
10
−4
10
−6
0
0
10
sample value − Sk
sample value − Sk
sample value − Sk
0
10
10
2
10
400 600 sample number − k
800
1000
10
0
−2
10
−4
10
−6
200
10
−6
200
400 600 sample number − k
800
1000
10
0
200
400 600 sample number − k
800
1000
37/51
MAP: “Family Picture” - n ≤ 2 states Name Poisson Erlang Renewal Process Hyper-exp. Renewal Process Interrupted Poisson Process MMPP Acyclic MAP(2) MAP(2)
D0 −λ −λ λ 0 −λ −λ1 0 0 −λ2 −λ1 λ1,2 λ −λ2 2,1 −λ1 λ1,2 λ −λ2 2,1 −λ1 λ1,2 0 −λ2 −λ1 λ1,2 λ2,1 −λ2
D1 λ 0 0 λ 0 pλ1 qλ1 pλ2 qλ2 ∗ λ1,1 λ∗1,2 0 0 ∗ λ1,1 0 0 λ∗2,2 ∗ λ1,1 0 λ∗ λ∗2,2 ∗2,1 λ1,1 λ∗1,2 λ∗2,1 λ∗2,2
Remarks: p + q = 1
38/51
MAP: “Family Picture” - Examples – n = 3 states Name −λ1 0 Hyper-Exp. Renewal 0 −λ1 Circulant MMPP(3) λ2,1 0 −λ1 λ2,1 MMPP(3) λ3,1 −λ1 0 Hyper-Exp. MAP 0 −λ1 λ2,1 MAP(3) λ3,1 Remarks: p + q + r = 1
D0 0 −λ2 0 0 −λ2 λ3,2 λ1,2 −λ2 λ3,2 0 −λ2 0 λ1,2 −λ2 λ3,2
0 0 −λ3 λ1,1 0 −λ3 λ1,3 λ2,3 −λ3 0 0 −λ3 λ1,3 λ2,3 −λ3
λ1 p λ1 p λ p ∗3 λ1,1 0 0 ∗ λ1,1 0 0 ∗ λ1,1 λ∗2,1 λ∗ 3,1 λ∗1,1 λ∗2,1 λ∗3,1
D1 λ1 q λ2 q λ3 q 0 λ∗2,2 0 0 λ∗2,2 0 λ∗1,2 λ∗2,2 λ∗3,2 λ∗1,2 λ∗2,2 λ∗3,2
λ1 r λ2 r λ3 r 0 0 λ∗3,3 0 0 λ∗3,3 λ∗1,3 λ∗2,3 λ∗3,3 λ∗1,3 λ∗2,3 λ∗3,3
39/51
MAP: Key Formulas for Inter-Arrival Times
• MAP representation: (D0 , D1 ) • Embedded chain: P = (−D0 )−1 D1 • Interval-stationary initial vector: α = αP Distribution of samples: • (Marginal) Distribution: F (t) = 1 − αe D0 t 1
• (Marginal) Density: f (t) = αe D0 t (−D0 )1 = αe D0 t D1 1
• (Marginal) Moments: E [X k ] = k!α(−D0 )−k 1 • Degrees of freedom of distribution: 2n − 1
40/51
MAP: Key Formulas for Inter-Arrival Times Sequence of samples: • Joint Density Function: Pr(X1 = x1 , X2 = x2 , . . . , Xq = xq ) = αe D0 x1 D1 e D0 x2 D1 · · · e D0 xq D1 1
• D1 → D1 Ph−1 for samples spaced by h lags • Joint Moments: E [X1k1 X2k2 · · · Xq q ] = k1 ! · · · kq !α(−D0 )−k1 P(−D0 )−k2 · · · (−D0 )−kq 1 k
• P → Ph for samples spaced by h lags • Analysis often limited to second-order properties. • Autocorrelation function: ρh =
α(−D0 )−1 Ph (−D0 )−1 1 − α(−D0 )−1 1 E [X1 X1+h ] − E [X ]2 = E [X 2 ] − E [X ]2 2α(−D0 )−2 1 − α(−D0 )−1 1 References: [TelH07]
41/51
4. INTER-ARRIVAL PROCESS FITTING
Includes joint work with Evgenia Smirni (IBM Research, William & Mary)
42/51
MAP: Approximate Inter-Arrival Process Fitting KPC readily generalizes to MAP random variables N (j) N (j) N • X = Jj Xj ∼ ((−1)J−1 Jj D0 , Jj D1 ) (j)
• X is a MAP if J − 1 D0 are diagonal (hyper-exp. moments) Moments and Joint Moments: Q E [X k ] • E [X k ] = k! Jj=1 k!j u QJ E [Xj,ik Xj,i+h ] u ] = k!u! • E [Xik Xi+h , ∀ lags i, h j=1 k!u! u v QJ E [Xj,ik Xj,i+h Xj,i+h+l ] u Xv • E [Xik Xi+h , ∀ lags i, h, l j=1 i+h+l ] = k!u!v ! k!u!v ! • ... References: [CasZS07]
43/51
MAP: KPC-Toolbox Autocorrelation decay rates and embedded chain: Q N • γi = Jj=1 γj,i , P = Jj=1 Pj Second-order properties follow recursively from those of X ⊗ Y : 1 + c 2 = 21 (1 + cX2 )(1 + cY2 ) 2 2 2 2 cX cY cX cY X Y Y ρk = ρk + ρk + ρX k ρk 2 2 c c c2 KPC-Toolbox – MAP Fitting: • Optimization-based second-order and third-order fitting • Xj ∼ MAP(2), known feasibility region for c 2 and ρk • Fitting of c 2 and ρk disjoint from first-order and third-order • Residual df spent to fit third-order moments E [Xi Xi+h Xi+h+l ] References: [CasZS08]
44/51
Example 4.1: Second-Order vs Third-Order Fitting • Second-Order: approximately match E [Xi Xi+h ] (equiv. ρh ) • Third-Order: approx. match E [Xi Xi+h ] and E [Xi Xi+h Xi+h+l ] LBL−PKT−5
LBL−PKT−5
Pr(queue ≥ x) [log]
Trace Autocorrelation Identical Fittings
−1
−2
k
ρ IAT autocorrelations [log]
0
Trace Simulation Second−Order Fit Third−Order Fit
0
−2
−4
−3 0
1
2 3 lag k [log]
4
0
2 x [log]
4
References: [AndN02],[CasZS07]
45/51
Example 4.2: KPC MAP Fitting Comparison of inter-arrival and counting process methods: • KPC: second-order and third-order fit of inter-arrival process • Superposition: fit Hurst coefficient for counts SWeb − MAP/D/1 − Utilization 80%
SWeb − Autocorrelation Fitting
0
0
−1 Pr(queue ≥ x) [log]
Autocorrelation ρk [log]
Web MAP(16) kpc −1
−2
−2 −3 −4 −5
−3 0
1
2 3 Lag k [log]
4
5
−6 0
Trace SUP KPC 2
4
6
x [log]
References: [AndN98],[CasZS07]
46/51
CONCLUSION
47/51
Summary
• PHs and MAPs are tractable models for characterizing empirical data using Markov chains • PHs and MAPs are compositional with other Markov models, e.g., queueing systems • PHs and MAPs closed under several operations (mixture, convolution, kronecker sum, kronecker product, ...) • Models with n ≤ 3 states are analytically tractable • Larger models can be dealt with using divide-and-conquer methods (KPC)
48/51
REFERENCES
49/51
References and Additional Readings 1/2 [AldS87] D. Aldous and L. Shepp. The least variable phase type distribution is Erlang. Comm. Stat.-Stochastic Models, 3(3):467–473, 1987. [AndN98] A. T. Andersen and B. F. Nielsen. A Markovian approach for modeling packet traffic with long-range dependence. IEEE JSAC, 16(5):719–732, 1998. [AndN02] A. T. Andersen and B. F. Nielsen. On the use of second-order descriptors to predict queueing behavior of MAPs. Naval Res. Log., 49(4):391–409, 2002. [BauBK10] Falko Bause, Peter Buchholz, Jan Kriege. ProFiDo - The Processes Fitting Toolkit Dortmund. QEST 2010: 87-96. [BuchK09] P. Buchholz, Jan Kriege. A Heuristic Approach for Fitting MAPs to Moments and Joint Moments. In Proc of QEST, 2009: 53-62, 2009. [BodHGT08] L. Bodrog and A. Heindl and G. Horv´ ath and M. Telek. A Markovian Canonical Form of Second-Order Matrix-Exponential Processes. European Journal of Operations Research, 190(457-477), 2008. [Bre78] J. Brewer. Kronecker products and matrix calculus in system theory. IEEE Trans. on Circuits and Systems, 25(9):772–781, 1978. [CasZS07] G. Casale, E. Zhang, and E. Smirni. Characterization of moments and autocorrelation in MAPs. ACM Perf. Eval. Rev., 35(1):27–29, 2007. [CasZS08] G. Casale, E. Zhang, and E. Smirni. KPC-toolbox: Simple yet effective trace fitting using markovian arrival processes. In Proc. of QEST, 83–92, 2008. [CasZS08b] G. Casale, E. Z. Zhang, and E. Smirni. Interarrival times characterization and fitting for Markovian traffic analysis. Technical Report WM-CS-2008-02, College of William and Mary, 2008.
50/51
References and Additional Readings 2/2 [HeiHG06] A. Heindl, G. Horvath, and K. Gross. Explicit Inverse Characterizations of Acyclic MAPs of Second Order. EPEW workshop, 108-122, 2006. [HeiML06] A. Heindl, K. Mitchell, and A. van de Liefvoort. Correlation bounds for second-order MAPs with application to queueing network decomposition. Performance Evaluation, 63(6):553–577, 2006. [HorT09] G. Horvath, M. Telek. On the canonical representation of phase type distributions. Performance Evaluation, 66(8):396–409, (2009). [MaiO92] R. S. Maier and C. A O’Cinneide. A closure characterisation of phase-type distributions. J. App. Probab., 29:92–103, 1992. [Neu89] M. F. Neuts. Structured Stochastic Matrices of M/G/1 Type and Their Applications. Marcel Dekker, New York, 1989. [Nie98] B. F. Nielsen. Note on the Markovian Arrival Process, Nov 1998. http://www2.imm.dtu.dk/courses/04441/map.pdf [OsoH06] T. Osogami, M. Harchol-Balter. Closed form solutions for mapping general distributions to quasi-minimal PH distributions. Performance Evaluation, 63(6): 524-552, 2006. [TelH03] M. Telek and A. Heindl. Matching moments for acyclic discrete and continuous phase-type distributions of second order. International Journal of Simulation, 3:47–57, 2003. [TelH07] M. Telek and G. Horv´ ath. A minimal representation of Markov arrival processes and a moments matching method. Performance Evaluation, 64(9-12):1153-1168, 2007.
51/51
