![TABLE I ESTIMATED TECHNOLOGY PARAMETERS AND 3σ VARIATIONS [1]](/figures/table-i-estimated-technology-parameters-and-3s-variations-1-3immfqkc.webp)
![Fig. 10. Probability density function of the log-normal distribution ∆x ∼ exp[N(0, 1/3)], i.e., the logarithm of ∆x is a normal distribution whose mean is 0 and standard deviation is 1/3.](/figures/fig-10-probability-density-function-of-the-log-normal-276qrl4s.webp)
16 IEEE TRANSACTIONS ON COMPUTER-AIDED DESIGN OF INTEGRATED CIRCUITS AND SYSTEMS, VOL. 26, NO. 1, JANUARY 2007
Asymptotic Probability Extraction for Nonnormal
Performance Distributions
Xin Li, Member, IEEE, Jiayong Le, Student Member, IEEE, Padmini Gopalakrishnan, Student Member, IEEE,
and Lawrence T. Pileggi, Fellow, IEEE
Abstract—While process variations are becoming more signifi-
cant with each new IC technology generation, they are often mod-
eled via linear regression models so that the resulting performance
variations can be captured via normal distributions. Nonlinear re-
sponse surface models (e.g., quadratic polynomials) can be utilized
to capture larger scale process variations; however, such models
result in nonnormal distributions for circuit performance. These
performance distributions are difficult to capture efficiently since
the distribution model is unknown. In this paper, an asymptotic-
probability-extraction (APEX) method for estimating the un-
known random distribution when using a nonlinear response
surface modeling is proposed. The APEX begins by efficiently
computing the high-order moments of the unknown distribution
and then applies moment matching to approximate the character-
istic function of the random distribution by an efficient rational
function. It is proven that such a moment-matching approach
is asymptotically convergent when applied to quadratic response
surface models. In addition, a number of novel algorithms and
methods, including binomial moment evaluation, PDF/CDF shift-
ing, nonlinear companding and reverse evaluation, are proposed
to improve the computation efficiency and/or approximation ac-
curacy. Several circuit examples from both digital and analog
applications demonstrate that APEX can provide better accuracy
than a Monte Carlo simulation with 10
4
samples and achieve
up to 10× more efficiency. The error, incurred by the popular
normal modeling assumption for several circuit examples designed
in standard IC technologies, is also shown.
Index Terms—Circuit performance, probability, process var-
iation, response surface modeling.
I. INTRODUCTION
A
S IC TECHNOLOGIES are scaled to the deep submi-
crometer region, process variations are becoming critical
and significantly impact the overall performance of a circuit.
Table I shows some typical process parameters and their 3σ
variations as technologies are scaled from 250 to 70 nm.
These large-scale variations introduce uncertainties in circuit
behavior, thereby making IC design increasingly difficult. Low
product yield or unnecessary overdesign cannot be avoided if
Manuscript received June 22, 2005; revised December 29, 2005 and March
29, 2006. This work was supported in part by the MARCO Focus Center for Cir-
cuit & System Solutions (C2S2, www.c2s2.org) under Contract 2003-CT-888.
This paper was presented in part at the IEEE/ACM International Conference
on Computer Aided Design (ICCAD) 2004. This paper was recommended by
Associate Editor J. R. Phillips.
The authors are with the Department of Electrical and Computer
Engineering, Carnegie Mellon University, Pittsburgh, PA 15213 USA
(e-mail: xinli@ece.cmu.edu; jiayongl@ece.cmu.edu; pgopalak@ece.cmu.edu;
pileggi@ece.cmu.edu).
Digital Object Identifier 10.1109/TCAD.2006.882593
TABLE I
E
STIMATED TECHNOLOGY PARAMETERS AND 3σ VARIATIONS [1]
process variations are not accurately modeled and analyzed
within the IC design flow.
During the past decade, various statistical analysis tech-
niques [1]–[7] have been proposed and utilized in many appli-
cations such as statistical timing analysis, mismatch analysis,
yield optimization, etc. The objective of these techniques is to
model the probability distribution of the circuit performance
under random process variations. Nassif [1] applies a linear
regression to approximate a given circuit performance f (e.g.,
delay, gain, etc.) as a function of the process variations (e.g.,
∆V
TH
, ∆T
OX
, etc.) and assumes that all random variations
are normally distributed. As such, the performance f is also
a normal distribution, since the linear combination of normally
distributed random variables is still a normal distribution [8].
The linear regression model is efficient and accurate when
process variations are sufficiently small. However, the large-
scale variations in deep submicrometer technologies, which
reach almost ±50% in Table I, suggest the need for higher
order regression models in order to guarantee high approxima-
tion accuracy [4]–[7]. Using a higher order response surface
model, however, brings about new challenges due to the non-
linear mapping between the process variations and the circuit
performance f. The distribution of f is no longer normal, unlike
the case of the linear model. The authors in [3]–[5] utilize the
Monte Carlo simulation to evaluate the probability distribution
of f, but this is computationally expensive. Note that reducing
the computational cost for this probability extraction is crucial,
especially when the extraction procedure is an inner loop within
an optimization flow.
In this paper, we propose a novel asymptotic-probability-
extraction (APEX) approach for estimating the unknown
random distribution using the nonlinear response surface
modeling. Given a circuit performance f (e.g., the delay of
a digital circuit path or the gain of an analog amplifier), the
response surface modeling approximates f as a polynomial
function of the process variations (e.g., ∆V
TH
, ∆T
OX
, etc.).
Since the process variations are modeled as random variables,
0278-0070/$25.00 © 2007 IEEE
LI et al.: ASYMPTOTIC PROBABILITY EXTRACTION FOR NONNORMAL PERFORMANCE DISTRIBUTIONS 17
Fig. 1. Overall flow of APEX.
the circuit performance f, which is a function of these random
variables, is also a random variable. The APEX applies a
moment matching to approximate the characteristic function
of f (i.e., the Fourier transform of the probability density
function (PDF) [8]) by a rational function H. We conceptually
consider H to be of the form of the transfer function of a linear
time-invariant (LTI) system, and the PDF and the cumulative
distribution function (CDF) of f are approximated by the
impulse response and the step response of the LTI system H,
respectively. Fig. 1 shows the overall flow of the APEX. In this
paper, we assume that the response surface model is already
available and it is provided to the APEX for estimating the
probability distribution of the circuit performance.
We prove that the moment-matching approach utilized in
APEX is asymptotically convergent when applied to quadratic
response surface models. In other words, given a quadratic re-
sponse surface model f that is a quadratic function of normally
distributed random variables, the PDF and the CDF of f can
be uniquely determined by its moments {m
k
,k =1, 2,...,K}
when K approaches infinity, i.e., K → +∞. The PDF and
the CDF extracted by the APEX can be used to characterize
and/or optimize the statistical performance of analog and digital
circuits under process variations.
APEX extends existing moment-matching methods via four
important new contributions which significantly reduce the
computational cost and improve the approximation accuracy
for this particular application. First, a key operation required
by the APEX is the computation of the high-order moments,
which is extremely expensive when using traditional techniques
such as the direct moment evaluation. In APEX, we propose a
binomial evaluation scheme to recursively compute the high-
order moments for a given quadratic response surface model.
The binominal-moment-evaluation scheme is derived from a
statistical independence theory and achieves a speedup of more
than 10
6
× over the direct moment-evaluation technique in our
tested examples.
Second, the APEX approximates the unknown PDF by the
impulse response of an LTI system. Directly applying such
an approximation to any circuit performance with negative
value is infeasible, since it results in an LTI system that is
noncausal. To overcome this difficulty, the APEX applies a
modified Chebyshev inequality for PDF/CDF shifting.
Third, in several practical applications, we observe that a
direct moment matching is not efficient when the unknown PDF
has a large negative skewness. In such cases, we propose to ap-
ply a nonlinear-companding scheme to automatically compress
and/or expand the PDF such that the transformed PDF almost
has a zero skewness and can be approximated accurately.
Finally, the best case performance (e.g., the 1% point on
CDF) and the worst case performance (e.g., the 99% point
on CDF) are two important metrics to be evaluated in many
practical applications. Direct moment matching cannot capture
the 1% point value accurately since the moment-matching
approximation is most accurate for low-frequency components
(corresponding to the final value of a CDF) and least accurate
for high-frequency components (corresponding to the initial
value of a CDF). To address this problem, a reverse-evaluation
technique is proposed in this paper to produce an accurate
estimation of the 1% point.
The remainder of this paper is organized as follows. In
Section II, we review the background on principal component
analysis (PCA), response surface modeling, and classical mo-
ment problem. We propose our APEX approach in Section III
and discuss several implementation issues, including binomial
moment evaluation, PDF/CDF shifting, nonlinear companding,
and reverse evaluation in Section IV. We extend the APEX
algorithm to handle nonnormal process variations in Section V
and discuss the applications of APEX in Section VI. The
efficacy of APEX is demonstrated by several circuit examples
in Section VII, followed by our conclusions in Section VIII.
II. B
ACKGROUND
A. PCA
The PCA [9] is a statistical method that finds a set of
independent factors to represent a set of correlated random
variables. Given N process parameters X =[x
1
,x
2
,...,x
N
]
T
,
the process variations ∆X = X − X
0
, where X
0
contains the
mean values of X, are often approximated as zero-mean normal
distributions, and the correlations of ∆X can be represented by
a symmetric positive semidefinite correlation matrix R.PCA
decomposes R as
R = V · Σ · V
T
(1)
where Σ = diag(λ
1
,λ
2
,...,λ
N
) contains the eigenvalues of
R and V =[V
1
,V
2
,...,V
N
] contains the corresponding eigen-
vectors that are orthonormal, i.e., V
T
V = I (I is the identity
matrix). Based on Σ and V , the PCA defines a set of new
random variables
∆Y =Σ
−0.5
· V
T
· ∆X. (2)
These new random variables ∆Y are called the principal com-
ponents or factors. It is easy to verify that ∆Y are independent
and satisfy the standard normal distribution N (0, 1) (i.e., zero
mean and unit standard deviation).
The essence of PCA can be interpreted as a coordinate
rotation of the space defined by the original random vari-
ables. In addition, if the magnitude of the eigenvalues {λ
i
}
decreases quickly, it is possible to use a small number of
random variables, i.e., a small subset of principal components,
to approximate the original N-dimensional space. More details
on PCA can be found in [9].
B. Response Surface Modeling
Given a circuit topology, the circuit performance (e.g.,
gain, delay, etc.) is a function of the design parameters
18 IEEE TRANSACTIONS ON COMPUTER-AIDED DESIGN OF INTEGRATED CIRCUITS AND SYSTEMS, VOL. 26, NO. 1, JANUARY 2007
(e.g., bias current, transistor sizes, etc.) and the process pa-
rameters (e.g., V
TH
, T
OX
, etc.). The design parameters are
optimized and fixed during the design phase; however, the
process parameters must be modeled as random variables to
account for any random variations. After the PCA shown in
(1) and (2), the process variations can be represented by N
independent random variables ∆Y that satisfy the standard nor-
mal distribution N(0, 1). Therefore, given a set of fixed design
parameters, the circuit performance f can be approximated by
a linear regression model of ∆Y [1]
f(∆Y )=B
T
· ∆Y + C (3)
where B ∈ R
N
stands for the linear coefficients and C ∈ R
is the constant term. The linear regression model in (3) is
accurate when process variations are small. However, such a
linear model is becoming increasingly inaccurate as process
variations become relatively large in nanoscale technologies.
In addition, the linear regression model in (3) always yields
the conclusion that the worst case performance appears at one
of the process corners. This conclusion might not be valid
for many circuit performances, especially those analog circuit
performances with strong nonlinearities, e.g., the offset voltage
of an analog amplifier. For these reasons, applying quadratic re-
sponse surface models might be required to provide a sufficient
modeling accuracy [3]–[5] as
f(∆Y )=∆Y
T
· A · ∆Y + B
T
· ∆Y + C (4)
where A ∈ R
N×N
denotes the quadratic coefficients, B ∈ R
N
represents the linear coefficients, and C ∈ R is the constant
term. Without loss of generality, we assume that A is symmetric
in this paper, since any asymmetric quadratic form X
T
AX
can be easily converted to an equivalent symmetric form 0.5 ·
X
T
(A + A
T
)X [10].
The model coefficients in (3) and (4) can be determined
by solving a set of overdetermined linear equations over a
number of sampling points [11]–[13]. These sampling points
are typically generated from SPICE simulations or measure-
ment results. Since generating each sampling point can be
quite expensive, many response surface modeling algorithms
[11]–[13] attempt to minimize the required number of sampling
points while maintaining good modeling quality. For this pur-
pose, a great number of algorithms for design of experiments
(DOE) have been proposed [14], [15]. Many of these DOE
algorithms (e.g., fractional factorial design, orthogonal array,
etc.) generate deterministic sampling points. These sampling
points do not represent the actual probability distributions of
process variations and, therefore, cannot be used for estimating
performance distributions. Even if a DOE algorithm (e.g., Latin
hypercube sampling [15]) generates random sampling points,
the number of these sampling points is typically too small
(e.g., 10–100) to accurately estimate PDF/CDF functions. For
this reason, after the response surface model is created, special
techniques are required to extract the probability distribution.
C. Classical Moment Problem
The classical moment problem was first proposed and studied
by T. Stieltjes in 1894. The Stieltjes moment problem is defined
as follows [16].
Definition 1 (Stieltjes Moment Problem): Given a sequence
of numbers {m
k
,k =1, 2,...}, find a nondecreasing function
F (f), where f ∈ [0, +∞), such that
m
k
=
+∞
0
f
k
dF (f). (5)
Note that if the symbol f in (5) represents a random vari-
able, the numbers {m
k
,k =1, 2,...} are its moments and the
function F (f) is its CDF. In the Stieltjes moment problem, the
random variable f is restricted to be nonnegative, i.e., its CDF
F (f) is defined only for f ≥ 0. By varying the interval in which
F (f) is valid, two further variations of the classical moment
problem can be defined [16].
Definition 2 (Hamburger Moment Problem): Given a se-
quence of numbers {m
k
,k =1, 2,...}, find a nondecreasing
function F (f), where f ∈ (−∞, +∞), such that
m
k
=
+∞
−∞
f
k
dF (f). (6)
Definition 3 (Hausdorff Moment Problem): Given a se-
quence of numbers {m
k
,k =1, 2,...}, find a nondecreasing
function F (f), where f ∈ [0, 1], such that
m
k
=
1
0
f
k
dF (f). (7)
As shown in (6) and (7), the Hamburger moment prob-
lem and the Hausdorff moment problem are defined in the
interval (−∞, +∞) and [0,1], respectively. In this paper, the
process variations are modeled as normal distributions, which
are unbounded and distributed over (−∞, +∞). Therefore, the
circuit performance approximated by the quadratic model (4)
is also unbounded. The probability extraction problem that we
aim to solve is the Hamburger moment problem in (6).
The classical moment problem has been widely studied by
mathematicians for over 100 years, focusing on the theoretical
aspects of the problem, e.g., the existence and uniqueness of
the solution. Details of these theoretical results can be found
in [16] or other recent publications, e.g., in [17]. However,
the practical applications of this moment problem, especially
the computation efficiency of solving the problem, have not
been sufficiently explored. In this paper, we develop the APEX
algorithm which aims to solve the moment problem efficiently,
i.e., to improve the approximation accuracy and reduce the
computational cost for practical applications.
III. APEX
Given the quadratic response surface model in (4), the ob-
jective of APEX is to estimate the PDF pdf(f) and the CDF
LI et al.: ASYMPTOTIC PROBABILITY EXTRACTION FOR NONNORMAL PERFORMANCE DISTRIBUTIONS 19
cdf(f) for the performance f.
1
Instead of running expensive
Monte Carlo simulations, the APEX tries to find an Mth
order LTI system H whose impulse response h(t) and step
response s(t) are the optimal approximations for pdf(f) and
cdf(f), respectively.
2
The optimal approximation is determined
by matching the first 2M moments between h(t) and pdf(f) for
an Mth order approximation. In this section, we first describe
the mathematical formulation of the APEX algorithm. Then,
we link the APEX to probability theory and explain why it
is efficient in approximating PDF/CDF functions. Finally, we
prove that the moment-matching method utilized in APEX is
asymptotically convergent when applied to quadratic response
surface models.
A. Mathematical Formulation
We define the time moments [18] for a given circuit perfor-
mance f whose PDF is pdf(f) as follows:
s
k
=
(−1)
k
k!
+∞
−∞
f
k
· pdf(f) · df . (8)
In (8), the definition of time moments is identical to the tradi-
tional definition of moments in probability theory except for the
scaling factor (−1)
k
/k!.
Similarly, the time moments can be defined for an LTI
system H [18]. Given an Mth order LTI system whose transfer
function and impulse response are
H(s)=
M
i=1
a
i
s − b
i
and h(t)=
M
i=1
a
i
e
b
i
t
, (if t ≥ 0)
0, (if t<0).
(9)
The time moments of H are defined as [18]
s
k
=
(−1)
k
k!
+∞
−∞
t
k
· h(t) · dt = −
M
i=1
a
i
b
k+1
i
. (10)
In (9), the poles {b
i
,i=1, 2,...,M} and residues {a
i
,i=
1, 2,...,M} are the 2M unknowns that need to be determined.
Matching the first 2M moments in (8) and (10) yields the
following 2M nonlinear equations:
−
a
1
b
1
+
a
2
b
2
+ ···+
a
M
b
M
= s
0
−
a
1
b
2
1
+
a
2
b
2
2
+ ···+
a
M
b
2
M
= s
1
.
.
.
.
.
.
−
a
1
b
2M
1
+
a
2
b
2M
2
+ ···+
a
M
b
2M
M
= s
2M−1
. (11)
1
In this paper, ∆X, ∆Y ,and∆Z represent the random variables for mod-
eling process variations, and f represents the circuit performance of interest.
The probability density function and the cumulative distribution function are
functions of f. Therefore, they are denoted as pdf(f) and cdf(f ), respectively.
2
The variable t in h(t) and s(t) corresponds to the variable f in pdf(f) and
cdf(f).
The nonlinear equations in (11) can be solved using the algo-
rithm proposed in [18], which first solves the poles {b
i
} and
then the residues {a
i
}. In what follows, we briefly describe this
two-step algorithm for solving (11).
1) Solving Poles: In order to solve the poles {b
i
} in (11),
Pillage and Rohrer [18] first formulate the following linear
equations:
−
s
0
s
1
··· s
M−1
s
1
s
2
··· s
M
.
.
.
.
.
.
.
.
.
.
.
.
s
M−1
s
M
··· s
2M−2
·
c
0
c
1
.
.
.
c
M−1
=
s
M
s
M+1
.
.
.
s
2M−1
.
(12)
After solving (12) for {c
i
,i=0, 1,...,M − 1}, the poles {b
i
}
in (11) are equal to the reciprocals of the roots of the following
characteristic polynomial:
c
0
+ c
1
b
−1
+ c
2
b
−1
+ ···+ c
M−1
b
−M+1
+ b
−M
=0. (13)
The detailed proof of (12) and (13) can be found in [18].
2) Solving Residues: After the poles {b
i
}are known, substi-
tute {b
i
} into (11) and the residues {a
i
} can be solved by using
the first M moments
−
b
−1
1
b
−1
2
··· b
−1
M
b
−2
1
b
−2
2
··· b
−2
M
.
.
.
.
.
.
.
.
.
.
.
.
b
−M
1
b
−M
2
··· b
−M
M
·
a
1
a
2
.
.
.
a
M
=
s
0
s
1
.
.
.
s
M−1
. (14)
The aforementioned algorithm assumes that the poles {b
i
}
are distinct. Otherwise, if repeated poles exist, the unknown
poles and residues must be solved using a more comprehensive
algorithm described in [18]. Once the poles {b
i
} and residues
{a
i
}are determined, the PDF pdf(f) is optimally approximated
by h(t) in (9) and the CDF cdf(f) is optimally approximated by
the step response:
s(t)=
t
0
h(τ)dτ =
M
i=1
a
i
b
i
· (e
b
i
t
− 1), (if t ≥ 0)
0, (if t<0).
(15)
It should be noted that many implementation issues must be
considered to make our proposed approach, APEX, feasible
and efficient. For example, the impulse response of a causal
LTI system is only nonzero for t ≥ 0, but a PDF in practical
applications can be nonzero for f ≤ 0. In Section IV, we will
propose several schemes to address these problems.
The aforementioned moment-matching method was previ-
ously applied to IC interconnect order reduction [18], [19],
and it is related to the Padé approximation in linear control
theory [20]. In the following section, we will explain why
such a moment-matching approach is efficient in approximating
PDF/CDF functions.
20 IEEE TRANSACTIONS ON COMPUTER-AIDED DESIGN OF INTEGRATED CIRCUITS AND SYSTEMS, VOL. 26, NO. 1, JANUARY 2007
Fig. 2. Characteristic functions of several typical distributions.
B. Connection to Probability Theory
In probability theory, given a random variable f whose PDF
is pdf(f), the characteristic function is defined as the Fourier
transform of pdf(f) [8]
Φ(ω)=
+∞
−∞
pdf(f) · e
jωf
· df =
+∞
−∞
pdf(f) ·
+∞
k=0
(jωf)
k
k!
· df .
(16)
Substituting (8) into (16) yields
Φ(ω)=
+∞
k=0
s
k
· (−jω)
k
. (17)
Equation (17) implies an important fact that the time moments
defined in (8) are related to the Taylor expansion of the char-
acteristic function at the expansion point ω =0. Matching the
first 2M moments in (11) is equivalent to matching the first 2M
Taylor expansion coefficients between the original characteris-
tic function Φ(ω) and the approximated rational function H(s).
To explain why the moment-matching approach is efficient,
we first need to show two important properties of the character-
istic function that are described in [8].
Property 1: A characteristic function has the maximal mag-
nitude at ω =0, i.e., |Φ(ω)|≤Φ(0) = 1.
Property 2: A characteristic function Φ(ω) → 0 when
ω →∞.
Fig. 2 shows the characteristic functions for several typical
random distributions. The above two properties imply an inter-
esting fact: Namely, given a random variable f, the magnitude
of its characteristic function decays as ω increases. Therefore,
the optimally approximated H(s) in (9) is a low-pass system.
It is well known that a Taylor expansion is accurate around the
expansion point. Since a low-pass system is mainly determined
by its behavior in the low-frequency range (around ω =0), it
can be accurately approximated by matching the first several
Taylor coefficients at ω =0, i.e., the moments. In addition, the
rational function form utilized in APEX is an efficient form to
approximate the transfer function H(s) of a low-pass system.
These conclusions have been verified in other applications
(e.g., IC interconnect order reduction [18], [19]), and they
provide the theoretical background to explain why the APEX
works well for PDF/CDF approximation, as will be demon-
strated by the numerical examples in Section VII.
C. Proof of Convergence
In the previous section, we have intuitively explained why
the moment-matching approach is efficient in approximating
PDF/CDF functions. However, there is a theoretical question
which might be raised: Given a random variable, can the
PDF/CDF functions always be uniquely determined by its
moments? In general, the answer is no. It has been observed
in mathematics that some probability distributions cannot be
uniquely determined by their moments. One example described
in [16] is the following PDF:
pdf(f)=
e
−0.5·[ln(f )]
2
√
2πf
·{1+a · sin [2π · ln(f)]}, (if f>0)
0, (if f ≤ 0)
(18)
where a ∈ [−1, 1]. It can be verified that all moments of the
PDF pdf(f) in (18) are independent of a, although varying a
changes pdf(f ) significantly [16]. It, in turn, implies that the
PDF in (18) cannot be uniquely determined by its moments.
However, there are special cases for which the moment prob-
lem is guaranteed to converge, i.e., the PDF/CDF functions are
uniquely determined by the moments. The following Carleman
theorem states one of those special cases and gives a sufficient
condition for the convergence of the moment problem.
Theorem 1 (Carleman [16]): A probability distribution on
the interval (−∞, +∞), i.e., the Hamburger moment problem
[see (6)], can be uniquely determined by its moments {m
k
,k =
1, 2,...} if
+∞
k=1
(m
2k
)
−1
2k
= ∞. (19)
Based on the Carleman theorem, we can prove that the moment-
matching approach utilized in APEX is asymptotically con-
vergent when applied to quadratic response surface models.
Namely, given a quadratic response surface model f that is a
quadratic function of the normally distributed random variables,
the probability distribution of f can be uniquely determined
by its moments {m
k
,k =1, 2,...,K} when K approaches
infinity, i.e., K → +∞. The asymptotic convergence of APEX
can be formally stated by the following theorem. The detailed
proof of Theorem 2 is given in Appendix.
Theorem 2: Given the quadratic response surface model f
in (4) where the random variables ∆Y are independent and
satisfy the standard normal distribution N(0, 1), the probability
distribution of f can be uniquely determined by its moments
{m
k
,k =1, 2,...}.
IV. I
MPLEMENTATIONS OF APEX
Our proposed APEX approach is made practically feasible
by applying several novel algorithms, including: 1) a binomial
scheme for high-order moment computation; 2) a modified