# Statistical Static Timing Analysis of VLSI as the Statistics of Correlated
Extremes

Dmytro Mishagli, Eugene Koskin, IEEE Member, and Elena Blokhina, Senior IEEE
Member
School of Electrical and Electronic Engineering, University College Dublin,
Ireland

###### Abstract

In this study, a path-based Statistical Static Timing Analysis (SSTA) is
formulated as a problem within the statistics of correlated extremes. For
extreme value statistics with correlations, a novel approach to studying such
systems, when the correlations are small, is developed. The approach considers
a system from the first principals, starting from the multivariable
characteristic function. Analytical solutions to the problem of weakly
correlated extremes are obtained in the form of corrections to the Gumbel
distribution. These solutions are compared with Monte Carlo simulations. The
applicability limits of the proposed solutions are studied. An algorithm to
estimate the covariance matrix of a timing graph is proposed.

###### Index Terms:

Timing Graph, Statistical Timing Analysis (SSTA), Delay, Very Large Scale
Integration (VLSI), Extreme Value Distribution

##  I Introduction

The timing verification of digital integrated circuits (ICs) is an absolutely
essential step: scaling down analogue and digital ICs poses very strict
restrictions on their operation, in particular, in terms of timing and delay.
Thus, timing analysis (TA) tools are used in order to verify and optimise
digital design before fabrication [1, 2, 3, 4, 5, 6, 7]. However, with the
reduction of feature size, the loss of fabricated ICs only increases, since
the effect of uncertainties on the performance of circuits increases
dramatically.

Process variations can be _global_ (also known as inter-die or die-to-die
variations) [8, 6, 9, 10, 11]. Global variations include:

  * •

Variations that appear due to changes in the temperature or the supply of
chemical compounds used during the fabrication processes (deposition, etching,
oxidation, and others.)

  * •

Wafer-to-wafer variations in the doping level of the substrate resulting in
substrate resistivity variations.

  * •

Die-to-die geometry variations due to misalignments of photo masks or
lithography.

In addition, there are _local_ variations (intra-die) that may affect
integrated circuits [5, 12]. Local variations may appear because of random
dopant fluctuations, interface roughness and random grains in the
polycrystalline structure of the gate dielectric and in the metal gates. Local
variations cause a change in individual circuit components.

Since the size of the technology node in the previous generations of IC was
large enough, global variations were seen as a dominant source of variation in
IC in the past [13]. Local variations are becoming more and more important now
because the feature size of the most recent technologies decreased so much so
it has become comparable with atomic-scale dimensions.

A combinational circuit can be described by an equivalent _timing graph_ ,
which can be denoted as $G=(V,E)$. In this notation, $V$ is the vertex set,
and $E$ is the edge set of this graph. The vertices of a timing graph can be
connected to two different types of nodes (edges), and they correspond to
different types of delays arising in a combinational circuit.

An example of a generic combinational logic circuit and its timing graph
obtained by converting the paths and gates into a set of vertices and edges is
shown in Fig. 1. In this timing graph, one can identify multiple paths from
the source to the sink of the circuit. Timing analysis deals with the analysis
of the timing graph and its paths. There is a simple transform that converts a
default timing graph into a one with the single source and single sink.

Returning to the verification of ICs, their timing graphs and timing analysis,
we recall that the nodes in a graph represent delays in a signal. The time
when a signal arrives at a gate (node) in a circuit can vary due to many
reasons. The input data may experience errors, the temperature and voltage
supplying a logic circuit may change, and there are manufacturing differences.
Hence, the main goal of _Static Timing Analysis_ is to verify that despite
these possible variations, all signals will arrive neither too early nor too
late, and hence proper circuit operation can be assured.

Figure 1: Example of a combinational logic circuit and the construction of its
timing graph with one source and sink.

In STA, the word _static_ implies that this timing analysis is carried out in
an input-independent manner. The purpose of this analysis is to find the
worst-case delay of a given circuit considering all possible input
combinations. It is known that the computational efficiency of STA is
excellent, and, hence, it is a common circuit verification tools, although it
has a number of limitations. In particular, STA cannot handle intra-die
variations and spatial correlation easily. If one wants to consider local
variations, STA needs many corners to handle all possible cases. In addition,
if there are significant random variations, STA predictions become too
pessimistic.

Known limitations of STA have led to the development of _SSTA_ , which
replaces deterministic timing of gates and interconnects with probability
distributions. The outcome of SSTA is a distribution of possible circuit
outcomes rather than a single outcome [14]. There are two main approaches with
the broad context of SSTA: _path-based_ and _block-based_ methods [4, 15].

A path-based approach looks into the total gate and interconnect delays on
specific paths. The statistical foundation of the method may look simple at a
glance. This approach requires that the paths of interest are identified prior
to starting the analysis. While it sounds like a straightforward technique, a
significant challenge may arise in the case of non-Gaussian correlated
distributions. Another potential issue is that some other paths, which are not
included in the original selection, may be relevant but not analysed. Hence,
path selection and correlation analysis are extremely important.

Closely related to the path-based SSTA are Monte Carlo (MC) methods. While MC
analysis of ICs is de facto golden standard for the designers, it has become
extremely expensive. Thus, there is ongoing research aiming to reducing the
simulation time while keeping MC’s accuracy. One of the problems of MC-based
methods is related to random samples generation [16, 17, 18, 19]. A quantile
extraction algorithm for an MC simulation was proposed in [20]. The usage of
order statistics [21] allowed one to decrease the number of samples required
for MC simulations. In [22, 23, 24], an approach to model variability based on
generalised lambda distributions is proposed.

A block-based approach looks into the arrival times and operation times for
each node (each logic gate). The advantage of this method is completeness as
it does not require any path selection. The main challenge is that the
statistical $\max$ and $\min$ operations, with and without the presence of
correlations, is needed, and it is recognised to be a hard technical problem
in Statistics. We refer the Reader to the review papers [25, 15, 26] for
further details on the block-based methods.

In this paper, the path-based SSTA problem is formulated as that of Extreme
Values Statistics. Then, the problem is relaxed to the case of weakly
correlated RVs. Analytical solutions to the problem of weakly correlated RVs
are obtained in the form of corrections to the Gumbel distribution. These
solutions are compared against MC numerical simulations. The applicability
limits of the proposed solutions are studies. The path-based approach to SSTA
formulated in this study requires the covariance matrix of a timing graph,
therefore, an algorithm to estimate the matrix is proposed. This research
continues and extends the results of the work [27], presented at the ISCAS
2019 (Sapporo, Japan).

The paper is organised as follows. In Section II, we outline the basics of
path-based SSTA, make necessary statistical preliminaries and present a
general statement of the problem. In Sections III and IV, we develop an
asymptotic theory for extreme value statistics, assuming that correlations are
weak. These sections somewhat repeat the work [27]. In Section V, the higher
terms of the expansion of CDF (PDF) of correlated extremes are considered.
Section VI discusses deviation from IID case. An algorithm for estimation of
the path covariance matrix is presented in Section VII. Finally, overall
discussion of the obtained results concludes this papers with Section VIII.

##  II Statement of the Problem in Path-Based SSTA

Consider a simple combinational logic circuit and its timing graph shown in
Fig. 1. In its most general formulation, the problem addressed within path-
based SSTA for such a circuit is to determine extreme delays arising in it.
This problem is equivalent to the problem of finding the maximum of a set of
$N$ _correlated_ RV:

$\zeta=\max(X_{1},X_{2},\ldots,X_{N}),$ (1)

where each RV component $X_{i}$ describes the accumulated (total) delay of the
relevant path of a circuit.

Since the number of nodes in timing graphs of modern VLSI systems is typically
of the order of $10^{9}$ and greater, see Ref. [9, 7], one can conclude that
the number of delay contributions in a given path is sufficiently large in
order to satisfy the requirements of the Central Limit Theorem. Thus, without
any loss of generality, we assume that all of the RVs $X_{i}$ that describe
the accumulated delay along the $i^{\text{th}}$ path in a graph $G$, are
distributed according to the Gaussian distribution (see Fig. 2).

Figure 2: Illustration of the statistics of extremes. On the top: a heatmap
with the probability density function of the RV $\zeta$, which is the maximum
of a set of $N$ standard Gaussian IID Random Variables,
$X_{i}\sim\mathcal{N}(0,1)$, against the number $N$. On the bottom: the man
value and the standard deviation for the same RV $zeta$. One can see that the
mean value increases (standard deviation decreases) with the increase of the
number $N$.

In the case of zero correlations between RVs, _i.e.,_ IID random variables,
the distribution of (1) is known and given by the Gumbel function. Let us
introduce specific notations. The CDF of the Gumbel distribution will be
denoted as follows:

$\Psi_{N}(z)=\exp\left(-e^{-\frac{z-\alpha}{\beta}}\right),$ (2)

where

$\alpha=\Phi^{-1}\left(1-\frac{1}{N}\right),\quad\beta=\frac{\sqrt{2\pi}}{N\cdot\varphi(\alpha)},$
(3)

and the subscript $N$ highlights the number of components in the RV vector.
The function $\Phi^{-1}(\cdot)$ is the inverse of the Gaussian CDF:

$\Phi(x)\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{1}{\sqrt{2\pi}}\int\limits_{-\infty}^{x}\varphi(x^{\prime})dx^{\prime}=\frac{1}{2}\left[1+\operatorname{erf}\left(\frac{x}{\sqrt{2}}\right)\right],$
(4)

where the standard definition of the error function is used:

$\operatorname{erf}(x)=\frac{2}{\sqrt{\pi}}\int\limits_{0}^{x}e^{-x^{\prime
2}}dx^{\prime}.$ (5)

The function $\varphi(\cdot)$ is referred here as the Gaussian kernel, and is
related to the normal distribution:

$\omega_{0}(x|\mu,\sigma)\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{1}{\sqrt{2\pi}\sigma}\varphi\left(\frac{x-\mu}{\sigma}\right),\quad\varphi(x)\stackrel{{\scriptstyle\text{def}}}{{=}}e^{-\frac{1}{2}x^{2}},$
(6)

The mean $\tilde{\mu}_{N}$ and variance $\tilde{\sigma}^{2}_{N}$ of the Gumbel
distribution are expressed in terms of the parameters $\alpha$ and $\beta$ as
follows:

$\tilde{\mu}_{N}=\alpha+\gamma\beta,\quad\tilde{\sigma}^{2}_{N}=\frac{\pi^{2}}{6}\beta^{2},$
(7)

where $\gamma\approx 0.5772$ is the Euler–Mascheroni constant. The PDF
$\psi_{N}(z)$ of the Gumbel distribution is given by the derivative of (2):

$\psi_{N}(z)=\frac{1}{\beta}\exp\left(-e^{-\frac{z-\alpha}{\beta}}-\frac{z-\alpha}{\beta}\right)$
(8)

In the case of the EVS of _correlated_ random variables (which is usually the
case for most of practical applications, including VLSI circuits) the
distribution of the extreme variable (1) is not known. However, there are
special cases that allow exact solutions, which give approximations to the
statistics of variable (1), namely strong and weak correlation cases.

For _strongly_ correlated RVs, there is no known universal approach to
undertake. However, there are particular cases and systems where solutions
were obtained by employing different methods. For instance:

  * •

Extreme statistics in random matrix theory [28, 29, 30]. The joint PDF of the
eigenvalues of a random Gaussian $N\times N$ matrix is known [31], and can be
interpreted as a Gibbs-Boltzmann weight for a Coulomb gas [32]. The challenge
was to study the fluctuations of the largest eigenvalue when
$N\rightarrow\infty$. It was shown that it is given by the so-called _Tracy-
Widom distribution_ 111In [33], this distribution was referred to as “the most
exciting recent result in mathematical physics”. [28, 29]. One should note
that the same distributions appear in quite different problems: in finance
[33], growth models [34], polymers [35, 36], mesoscopic fluctuations in
quantum dots [37], and many others.

  * •

Extreme value statistics in the presence of a global constraint. Such a
constraint (_e.g._ , it can be a conservation law) results in strong
correlations between RVs. An example that allowed a solution is known as the
zero range process [38], which is a 1D lattice of $N$ interacting agents with
periodic boundary conditions. The global constraint (and, therefore,
correlations) is realised mathematically via $\delta$-functions.

The detailed discussion of these and other cases can be found in the review
papers [39, 40], as well as in the literature cited there.

Figure 3: Illustration of a set of weakly correlated RVs studied within the
Renormalisation Group approach [41, 42]

The system of _weakly_ correlated RVs is traditionally studied via the
Renormalisation Group (RG) approach [41, 42]. The idea is illustrated in Fig.
3 and is as follows. For a set of $N$ correlated RVs, consider _blocks_ in
such a way, so that the RVs become non-correlated between the blocks. This is
possible to achieve if the correlation between the variables decays
exponentially, and is negligible after some characteristic correlation length
$\Delta$:

$\langle X_{i}X_{j}\rangle-\langle X_{i}\rangle\langle X_{j}\rangle\sim
e^{-\frac{|i-j|}{\Delta}}.$ (9)

Here, the angled brackets $\langle\cdot\rangle$ denote the averaging
operation. Thus, one obtains a set of new _uncorrelated_ RVs,
$y_{1},y_{2},\ldots,y_{N^{\prime}}$, where $N^{\prime}=N/\Delta$. Each $y_{i}$
represents the maximum in the corresponding block. The problem, however, lies
in determining this maximum within each block, as RVs within blocks remain
correlated strongly. Therefore, the methods of estimating
$y_{1},y_{2},\ldots,y_{N^{\prime}}$ should be applied depending on the system
under consideration, and the problem reduces to the strongly correlated cases
considered above.

In this study, for the purposes of SSTA, we consider the problem of finding
the distribution of the maximum of weakly correlated Gaussian RVs. The problem
is stated in the most general form: starting from the _first principles_ , we
aim to find the expression for the PDF of correlated Gaussian RVs and the
maximum of a set of RVs. This allows us to obtain corresponding corrections to
the Gumbel law (2), and this constitutes the original research contribution of
this paper.

##  III Weak Correlation Approximation

In this section, we start by deriving an expansion for the multivariate PDF
$\omega(\mathbf{r})$ of $N$ weakly-correlated Gaussian RVs $X_{i}$. We will
use this expansion in later Sections to obtain the main results of this work.
Let us use the notation of the random vector
$\mathbf{r}=(x_{1},x_{2},\ldots,x_{N})$ whose components are the realisations
of random variables $X_{i}$. In order to obtain the expansion for
$\omega(\mathbf{r})$, we will use the concept of characteristic functions.
Since we need to distinguish between many types of functions and also between
correlated and uncorrelated random variables, the following symbols will be
used. The functions $\omega(\mathbf{r})$ and $\chi(\mathbf{k})$ will denote
the PDF and characteristic function of _correlated_ Gaussian RVs respectively,
while the _uncorrelated_ case will be denoted using the same symbols but with
the subscript ‘0’.

The PDF $\omega(\mathbf{r})$, in the general case, is presented as follows:

$\omega(\mathbf{r})=\frac{1}{(2\pi)^{N}}\int\limits_{E_{N}}\chi(\mathbf{k})e^{-i\mathbf{k}\mathbf{r}}d\mathbf{k}\,,$
(10)

with the characteristic function

$\chi(\mathbf{k})=\exp\left(i\mu_{i}k_{i}-\frac{1}{2}\Sigma_{ij}k_{i}k_{j}\right).$
(11)

Here, $\Sigma_{ij}$ is the covariance matrix such that

$\Sigma_{ij}=\delta_{ij}\sigma_{i}^{2}+\varepsilon_{ij},\quad\varepsilon_{ii}=0.$
(12)

In the above equations, $\mu_{i}$ and $\sigma_{i}$ are the mean and standard
deviation of the RV components, which are assumed to be known. The latter
allows us to factorise the exponent in (11):

$\chi(\mathbf{k})=\chi_{0}(\mathbf{k})\cdot\exp\left(-\frac{1}{2}\varepsilon_{ij}k_{i}k_{j}\right),$
(13)

where
$\chi_{0}(\mathbf{k})=\exp\left(i\mu_{i}k_{i}-\frac{1}{2}\sigma_{i}^{2}k_{i}^{2}\right)$
is the characteristic function of $N$ uncorrelated normal RVs.

Let us assume that the non-diagonal entries of the covariance matrix are small
by absolute value, and, therefore, the characteristic function can be
presented through an expansion with the coefficients $\varepsilon_{ij}$ acting
as perturbations:

$\chi(\mathbf{k})\bigg{|}_{\varepsilon_{ij}\rightarrow
0}=\chi_{0}(\mathbf{k})+\varepsilon_{ij}\cdot\left.\frac{\partial\chi(\mathbf{k})}{\partial\varepsilon_{ij}}\right|_{\varepsilon_{ij}=0}+\dots$
(14)

We can express now the PDF $\omega(\mathbf{r})$ via the non-correlated term
$\omega_{0}$ with the following remarkable property of the multivariate normal
distribution:

$\left.\frac{\partial\chi}{\partial\varepsilon_{ij}}\right|_{\varepsilon_{ij}=0}=\frac{1}{2}\frac{\partial^{2}\chi_{0}}{\partial\mu_{i}\partial\mu_{j}}\,.$

Indeed, limiting ourselves to the first order of smallness and substituting
(14) into (10), we obtain the expansion for the PDF of normal correlated
random variables in the weak correlation case:

$\omega(\mathbf{r})\simeq\omega_{0}(\mathbf{r})+\frac{\varepsilon_{ij}}{2}\frac{\partial^{2}\omega_{0}(\mathbf{r})}{\partial\mu_{i}\partial\mu_{j}}\,.$
(15)

At this stage, we can introduce the explicit form for the PDF
$\omega_{0}(\mathbf{r})$. According to the note made at the beginning of this
Chapter, we assume that the CLT holds. Hence, the RVs in our case are seen as
correlated Gaussian ones. We can write $\omega_{0}(\mathbf{r})$ as the PDF of
$N$ Gaussian RVs with the mean $\mu_{k}$ and variance $\sigma_{k}$:

$\omega_{0}(\mathbf{r})\overset{\text{def}}{=}\frac{1}{(2\pi)^{N/2}}\prod\limits_{k=1}^{n}\frac{1}{\sigma_{k}}\varphi\left(\frac{x_{k}-\mu_{k}}{\sigma_{k}}\right).$
(16)

In the case of only one RV, we simply have (6).

In principle, the second derivative of $\omega_{0}(\mathbf{r})$ can be written
explicitly:

$\frac{\partial^{2}\omega_{0}(\mathbf{r})}{\partial\mu_{i}\partial\mu_{j}}=\frac{\omega_{0}(\mathbf{r})}{\sigma_{i}^{2}}\left[\frac{(x_{i}-\mu_{i})(x_{j}-\mu_{j})}{\sigma_{j}^{2}}-\delta_{ij}\right],$

where $\delta_{ij}$ is the Kronecker delta symbol. However, this step is not
necessary for the present discussion. Instead, we will proceed to
differentiation by components:

$\omega(\mathbf{r})\simeq\omega_{0}(\mathbf{r})+\frac{\varepsilon_{ij}}{2}\frac{\partial^{2}\omega_{0}(\mathbf{r})}{\partial
x_{i}\partial x_{j}}.$ (17)

It is easy to see that expressions (15) and (17) are identical. At the same
time, representation (17) allows one to treat the mean values, $\mu_{i}$ and
$\mu_{j}$, as parameters, which simplifies the analysis below. In the next
Section, we shall consider the CDF $\mathfrak{F}(x)$, assuming that the RVs
$X_{i}$ have the PDF in the form of (17).

##  IV Statistics of Weakly Correlated Extreme Values: First Order
Approximation

###  IV-A Derivation of the Cumulative Distribution Function in the First
Order Approximation

To find the CDF $\mathfrak{F}(z)$, we recall that the probability $\mathbb{P}$
that all of the RVs $X_{1},X_{2},\ldots,X_{N}$ are less than some given number
$z$, $\mathbb{P}[\zeta<z]$, is

$\mathbb{P}[\zeta<z]\equiv\mathfrak{F}(z)=\int\limits_{-\infty}^{z}\ldots\int\limits_{-\infty}^{z}\omega(\mathbf{r})d\mathbf{r},$
(18)

where $d\mathbf{r}=dx_{1}dx_{2}\ldots dx_{N}$ is the elementary volume in the
$N$-dimensional space of random variables $X_{i}$.

We have not imposed any restrictions on the mean values and standard
deviations, $\mu_{i}$ and $\sigma_{i}$, while deriving expansion (17). For the
sake of simplicity, we assume that $\mu_{i}$ and $\sigma_{i}$ are the same for
all $X_{i}$, _i.e.,_ we consider a set of IID random variables. In such a
case, RVs can be re-scaled to be standard normal ones,
$X_{i}\sim\mathcal{N}(0,1)$, using the transform:

$X_{i}\rightarrow\frac{X_{i}-\mu}{\sigma}.$

The effect of a deviations from the IID case is discussed below in the next
Section.

Substituting expansion (17) into (18) and letting $N\rightarrow\infty$, we get
two terms. The first term leads to the CDF of the Gumbel distribution:

$\int\limits_{-\infty}^{z}\ldots\int\limits_{-\infty}^{z}\omega_{0}(\mathbf{r})d\mathbf{r}=\Phi(z)^{N}\underset{N\rightarrow\infty}{=}\Psi_{N}(z),$

while the second one involves the integral of the form

$I=\int\limits_{-\infty}^{z}\ldots\int\limits_{-\infty}^{z}\frac{\partial^{2}\omega_{0}(\mathbf{r})}{\partial
x_{i}\partial x_{j}}d\mathbf{r}.$

The latter gives

$\displaystyle I$
$\displaystyle=\frac{1}{\sqrt{2\pi}}\int\limits_{-\infty}^{z}\frac{\partial}{\partial
x_{i}}e^{-\frac{x_{i}^{2}}{2}}dx_{i}\cdot\frac{1}{\sqrt{2\pi}}\int\limits_{-\infty}^{z}\frac{\partial}{\partial
x_{j}}e^{-\frac{x_{j}^{2}}{2}}dx_{j}$
$\displaystyle\times\frac{1}{(2\pi)^{(N-2)/2}}\int\limits_{-\infty}^{z}\ldots\int\limits_{-\infty}^{z}\exp\left(-\frac{1}{2}\sum\limits_{k\neq
i,j}x_{k}^{2}\right)\underbrace{dx_{1}\ldots dx_{N}}_{\text{except
}dx_{i}dx_{j}}$
$\displaystyle\underset{N\rightarrow\infty}{=}\frac{\varphi(z)^{2}}{2\pi}\Psi_{N-2}(z).$

Finally, bringing all terms together and taking into account that
$\Psi_{N-2}(x)\approx\Psi_{N}(x)$, we can write down the asymptotic expression
for the CDF of the maximum of $N$ weakly–correlated Gaussian random variables
in the first order approximation:

$\mathfrak{F}_{N}(z)=\Psi_{N}(z)\left[1+\frac{\varphi(z)^{2}}{4\pi}\sum_{\begin{subarray}{c}i,j\\\
i\neq j\end{subarray}}\varepsilon_{ij}+\ldots\right].$ (19)

The PDF of the extended Gumbel law (19) can be calculated explicitly. The
generalised formula for the PDF is given in the next Section by expression
(V-A) since it serves both cases, the first order and the second order
approximations.

Note that this formula is derived assuming the diagonal entries
$\varepsilon_{ii}=0$, and is valid for small correlations:

$|\varepsilon_{ij}|\ll\delta_{ij}\sigma_{i}^{2}.$ (20)

Formula (19) represents an extension of the Gumbel distribution and
constitutes the original result of this Section. Modelling and simulation of
this case are presented in the next Section.

###  IV-B Modelling and Comparison with Monte Carlo Simulations Using the
First Order Approximation

To understand the result and verify its validity, we will compare it with
Monte Carlo simulations. This type of simulations is the reference method to
verify a VLSI design [4]. For the sake of illustration of the obtained
results, we choose the correlation between the RVs $X_{i}$ as follows:

$\langle X_{i}X_{j}\rangle=\sigma^{2}\rho^{|i-j|},\quad\langle
Y_{i}Y_{j}\rangle=\delta_{ij},\quad\langle X_{i}Y_{j}\rangle=0,$ (21)

where $X_{i}\sim\mathcal{N}(0,\sigma)$ and $Y_{i}\sim\mathcal{N}(0,1)$. In
this case, random samples can be generated via the simple rule:

$X_{i+1}=\rho X_{i}+\sigma\sqrt{1-\rho^{2}}Y_{i},$ (22)

where $\rho$ is the correlation coefficient. A set of $N$ samples has been
generated according to rule (22) with the starting point $X_{0}$ chosen
randomly from $\sim\mathcal{N}(0,\sigma)$, then $\max(X_{i},\ldots,X_{N})$ was
computed. The procedure was repeated $10^{4}$ times for each numerical
experiment.

Figure 4: (These results were first presented in [27], at the ISCAS 2019, held
in Sapporo, Japan.) Illustration of the proposed first-order approximation
compared against the standard Gumbel distribution and Monte Carlo simulations.
The upper row shows the PDF of the weakly correlated Gaussian RVs (21) obtain
through the proposed approach and compared with the Gumbel distribution and MC
simulations. Two different examples of correlations $\rho=0.35$ and
$\rho=0.75$ are shown for the number of random variables $N=100$. The lower
graph shows the behaviour of the mean as a function of the correlation
coefficient $\rho$ for $N=200$.

For a large number of RVs ($N\gtrsim 10^{2}$), the dependence of the PDF on
the correlation coefficient $\rho$ is studied. As one can see from Fig. 4, for
relatively small values of $\rho$, the deviation of both the perturbed and
uncorrelated Gumbel distributions, (19) and (2) respectively, from MC
simulations is negligible (error in the mean value is less than 2$\%$). This
implies that small correlations can be ignored without any noticeable loss of
accuracy. At the same time, the comparison in the strong correlation regime,
$\rho>0.5$, shows that approximation (19), even though it should not be
applied for such $\rho$, is able to describe the correct trend.

Since this approach is developed with delay propagation in a logic circuit in
mind, it would be particularly interesting to discuss the behaviour of the
mean value as a function of the correlation coefficient $\rho$, as also shown
in Fig. 4. One can see that the uncorrelated Gumbel distribution (2) leads to
significant deviations from the numerical simulations for $\rho>0.5$. Thus,
the expression for the mean (7) can be interpreted as the worst case delay. At
the same time, weakly-correlated formula (19) gives a result between that and
the simulation dots. The deviation of (19) from the MC simulations may be
reduced by taking into account further terms in the expansion (14), as shown
in the next Section.

It should be pointed out that the mean is decreasing with the increase of the
correlation coefficient, which is seen both from the simulation and the theory
in Fig. 4. Recalling the meaning of the RV $\zeta$ in SSTA, we see that the
total delay in a logic circuit is inversely proportional to the value of the
correlation coefficient $\rho$.

##  V Statistics of Weakly Correlated Extreme Values: Higher Terms

###  V-A Derivation of the Cumulative Distribution Function in the Second
Order Approximation

In this Section, we will try to improve the previous result by taking into
account the second order term in the expansion (14):

$\displaystyle\chi(\mathbf{k})\bigg{|}_{\varepsilon_{ij}\rightarrow
0}=\chi_{0}(\mathbf{k})$
$\displaystyle+\varepsilon_{ij}\cdot\left.\frac{\partial\chi(\mathbf{k})}{\partial\varepsilon_{ij}}\right|_{\varepsilon_{ij}=0}$
$\displaystyle+\frac{1}{2}\varepsilon_{ij}\varepsilon_{kl}\cdot\left.\frac{\partial^{2}\chi(\mathbf{k})}{\partial\varepsilon_{ij}\partial\varepsilon_{kl}}\right|_{\varepsilon_{ijkl}=0}+\dots$
(23)

In the similar manner as was done for the first derivative, one can show that:

$\left.\frac{\partial^{2}\chi}{\partial\varepsilon_{ij}\partial\varepsilon_{kl}}\right|_{\varepsilon_{ijkl}=0}=\frac{1}{4}\frac{\partial^{2}\chi_{0}}{\partial\mu_{i}\partial\mu_{j}\partial\mu_{k}\partial\mu_{l}},$

and present the expansion of the PDF of correlated Gaussian RVs in the second
order approximation as follows:

$\omega(\mathbf{r})\simeq\omega_{0}(\mathbf{r})+\frac{\varepsilon_{ij}}{2}\frac{\partial^{2}\omega_{0}(\mathbf{r})}{\partial
x_{i}\partial
x_{j}}+\frac{\varepsilon_{ij}\varepsilon_{kl}}{8}\frac{\partial^{4}\omega_{0}(\mathbf{r})}{\partial
x_{i}\partial x_{j}\partial x_{k}\partial x_{l}}.$ (24)

Figure 5: Illustration of the proposed second-order approximation compared
against the first-order approximation, standard Gumbel distribution and Monte
Carlo simulations. The PDF is shown for four different correlation coefficient
$\rho=0.35$, $0.5$, $0.65$ and $0.75$, demonstrating the transition from
weakly-correlated to the strongly-correlated case. The number of variables is
$N=100$.

Now, substituting (24) into (18) and letting $N\rightarrow\infty$, we get the
term that involves integral of the following form:

$I=\int\limits_{-\infty}^{z}\ldots\int\limits_{-\infty}^{z}\frac{\partial^{4}\omega_{0}(\mathbf{r})}{\partial
x_{i}\partial x_{j}\partial x_{k}\partial
x_{l}}\underset{N\rightarrow\infty}{=}\frac{\varphi(z)^{4}}{(2\pi)^{2}}\Psi_{N-4}(z).$

Therefore, the asymptotic expression for the CDF of the maximum of $N$ weakly-
correlated Gaussian random variables in the second order approximation has the
form:

$\displaystyle\mathfrak{F}_{N}(z)=\Psi_{N}(z)\bigg{[}1$
$\displaystyle+\frac{\varphi(z)^{2}}{4\pi}\sum_{\begin{subarray}{c}i,j\\\
i\neq j\end{subarray}}\varepsilon_{ij}$
$\displaystyle+\frac{\varphi(z)^{4}}{32\pi^{2}}\sum_{\begin{subarray}{c}i,j,k,l\\\
i\neq j,k\neq
l\end{subarray}}\varepsilon_{ij}\varepsilon_{kl}+\ldots\bigg{]}.$ (25)

The corresponding expansion for the PDF $\mathfrak{f}_{N}(z)$ is obtained by
taking the derivative of (V-A). We have

$\displaystyle\mathfrak{f}_{N}(z)=\psi_{N}(z)$
$\displaystyle+\frac{\varphi(z)^{2}}{4\pi}a_{N}(z)\sum_{\begin{subarray}{c}i,j\\\
i\neq j\end{subarray}}\varepsilon_{ij}$
$\displaystyle+\frac{\varphi(z)^{4}}{32\pi^{2}}b_{N}(z)\sum_{\begin{subarray}{c}i,j,k,l\\\
i\neq j,k\neq l\end{subarray}}\varepsilon_{ij}\varepsilon_{kl}+\ldots,$ (26)

where the coefficients $a_{N}(z)$ and $b_{N}(z)$ read

$\displaystyle a_{N}(z)$ $\displaystyle=\psi_{N}(z)-2z\varphi(z)\Psi_{n}(z),$
(27) $\displaystyle b_{N}(z)$
$\displaystyle=\psi_{N}(z)-12z\varphi(z)^{3}\Psi_{n}(z).$

We note that in the PDF (V-A)-(27), the first perturbation term $a_{N}$
corresponds to the first order approximation considered in Section IV-A. The
second perturbation term $b_{N}$ appears due to the second order approximation
considered in this Section.

One can continue the expansion (V-A), and, letting $N\rightarrow\infty$, we
have

$\displaystyle\mathfrak{F}_{N}(z)$
$\displaystyle=\Psi_{N}(z)\sum_{n=0}^{\infty}\frac{1}{n!}\left(\frac{\varphi(z)^{2}S}{4\pi}\right)^{n}$
$\displaystyle=\Psi_{N}(z)\exp\left(\frac{\varphi(z)^{2}S}{4\pi}\right),$ (28)

where $S=\sum_{\begin{subarray}{c}i,j\\\ i\neq
j\end{subarray}}\varepsilon_{ij}$. Correspondingly, the PDF is given by the
derivative of (V-A).

Let us now examine the validity of the results obtained in this section and
compare them with the first order approximation (19).

Figure 6: Illustration of the proposed second-order approximation compared
against the first-order approximation, the standard Gumbel distribution and
Monte Carlo simulations. The mean value as a function of the correlation
coefficient $\rho$ is presented.

###  V-B Modelling and Comparison with Monte Carlo Simulations Using the
Second Order Approximation

The verification of the second order approximation has been carried out using
the same methodology outlined in Section IV-B. Modelling and comparison with
numerical simulations are summarised in Fig. 5. This figure illustrates the
proposed second-order approximation in action compared against the first-order
approximation, standard Gumbel distribution and Monte Carlo simulations. The
PDF of these four cases are shown for four different correlation coefficient
$\rho=0.35$, $0.5$, $0.65$ and $0.75$, demonstrating the transition from
weakly-correlated to strongly-correlated case. In all presented cases, one
sees a very good agreement between the PDF obtained from the second-order
approximation and the MC simulation.

The investigation of the dependence of the mean value on the correlation
coefficient $\rho$ is presented in Figure 6. The second order approximation
seems to be particularly effective for the intermediate number $N$ of the
components $X_{i}$. In all cases, both the first order and second order
approximations can capture the correct trend even in the presence of strong
correlations.

Figure 7: The obtained extensions to the Gumbel law, namely, the
$1^{\text{st}}$ and $2^{\text{nd}}$ approximations, and the complete formula
(V-A) for a case of $\rho=0.825$, where the weak correlation assumption is
clearly violated. Left column shows the case with $N=100$, and one can see how
the proposed corrections break down. Right column shows the case with $N=250$,
and the higher number of terms in the RV (1) leads to better performance, as
role of pure Gumbel dominates. The insets show $\log(-\log(\cdot))$ of the
CDF.

However, the obtained asymptotic results can quickly become inapplicable for
higher correlation coefficients. Hence Fig. 7 shows Gumbel CDF (2) and the
corrections (19), (V-A), and (V-A), as well as corresponding PDFs, for the
correlation coefficient $\rho=0.825$. There are two cases shown: relatively
small number of variables $N=100$ (left column), and the higher number of
$N=250$ (right column). For the first case, the asymptotic solutions blow
(compare with Fig. 5, where the PDFs are shown for $\rho\leq 0.75$ for the
same number of variables $N=100$). For the second case ($N=250$), the
corrections give good improvement of the Gumbel law. Such a behaviour is
expected, as the corrections are derived assuming $N\rightarrow\infty$.
Therefore, the applicability of the results to the realistic graphs with
smaller number of paths should be studied in a detail.

##  VI The Case of Non-Identical Independent Random Variables

The major approximation we made in the previous Sections on random variables
was to assume that they are identically distributed. While this assumption
simplified the analysis, in reality, one does not expect that the RVs will be
identical. If this approach is used to describe path delays, it not
necessarily the case that all nodes in logic circuits and delays accumulated
in sub-paths will have the same mean and variance. However, we highlight again
that due to the CLT, we still expect that the RVs will be Gaussian.

Figure 8: Effect of non-identical mean values and variances (standard
deviations) on the statistics of extreme values in the path-based approach.
The graphs shown on the left present the dependence of the mean and standard
deviation on the number $N$ of the random components $X_{i}$ in $\zeta$ when
one varies the strength of deviations $\delta\mu$ in the mean value. The
graphs on the right show their counterparts in the case when one varies the
strength of deviations $\delta\sigma$ in the standard deviation.

To investigate the effect of non-identical mean values and variances in the
random components $X_{i}$ comprising the variable $\zeta$, we consider small
deviations from IID. Namely, the mean values are assumed to be
$\mu_{i}=\mu+\xi\cdot\delta\mu$ while the variances are
$\sigma_{i}=\sigma+\xi\cdot\delta\sigma$. Here, the variable $\xi$ is another
RV drawn from the uniform distribution to help generate random deviations,
_i.e._ , $\xi\sim\mathcal{U}(-1,1)$. The quantities $\delta\mu$ and
$\delta\sigma$ give the _strength_ of deviations from the nominal $\mu$ and
$\sigma$. For simplicity, we assume that these non-identical RVs are
_uncorrelated_ or independent.

The effect of non-identical mean values and variances (standard deviations) on
the statistics of extreme values in the path-based approach is summarised in
Fig. 8. The graphs shown on the left present the dependence of the mean and
standard deviation on the number $N$ of the random components $X_{i}$ in the
variable $\zeta$ when one varies the strength of deviations $\delta\mu$ in the
mean value. The graphs on the right show their counterparts in the case when
one varies the strength of deviations $\delta\sigma$ in the standard
deviation. As a reference, the mean and standard deviation (2) calculated from
the Gumbel distribution of IID variables is shown by the black dashed line.
The only effect one can see in this case consists in the scaling of the
curves, keeping the overall shape of the Gumbel distribution. We conclude
that, if required, non-identical RVs can be accommodated in the approach
proposed in this paper by proper renormalisation of the relevant quantities.

##  VII Path Covariance Matrix

The final question we would like to cover in this Chapter is related to the
covariance matrix. Indeed, inspecting the example of the timing graph shown in
Fig. 1, we note that many paths starting at the Source and ending at the Sink
of the graph share common nodes. This is exactly the reason why correlations
appear in the random components $X_{i}$ in the path-based approach. Our task
is to propose a method to estimate the covariance matrix of a timing graph.

We will assume that the delays $\tau_{ij}$ of the nodes in a timing graph are
given by IID RVs drawn from a normal distribution, _i.e._ , $\tau_{ij}\sim
N(\mu,\sigma^{2})$, where $\mu$ is the mean value of a node delay and $\sigma$
is its standard deviation. In this case, the delay between the $i^{\text{th}}$
and $j^{\text{th}}$ nodes along a path can be written as

$\tau_{ij}=\mu+\sigma\xi_{ij},\quad\xi_{ij}\sim\mathcal{N}(0,1).$ (29)

Here, $\xi$ is an auxiliary random number drawn from the standard normal
distribution (with a zero mean and a unit variance) to help generate a set of
random delays. Thus,

$\langle\xi_{ij}\rangle=0,\quad\langle\xi_{ij}^{2}\rangle=1,\quad\langle\xi_{ij}\xi_{kl}\rangle=\delta_{ik}\delta_{jl}.$
(30)

Figure 9: An example of a directed acyclic graph illustrating shared nodes.

Such a representation allows one to estimate the covariance matrix due to
shared paths. For the sake of illustration, consider a simple directed acyclic
timing graph shown in Fig. 9. It is easy to see that there are four paths with
the accumulated delays as follows:

path 1: $\displaystyle\tau_{12}+\tau_{24}+\tau_{46}+\tau_{67}$
$\displaystyle=4\mu+\sqrt{4}\sigma\epsilon_{1}$ path 2:
$\displaystyle\tau_{12}+\tau_{24}+\tau_{43}+\tau_{35}+\tau_{56}+\tau_{67}$
$\displaystyle=6\mu+\sqrt{6}\sigma\epsilon_{2}$ path 3:
$\displaystyle\tau_{12}+\tau_{24}+\tau_{45}+\tau_{56}+\tau_{67}$
$\displaystyle=5\mu+\sqrt{5}\sigma\epsilon_{3}$ path 4:
$\displaystyle\tau_{13}+\tau_{35}+\tau_{56}+\tau_{67}$
$\displaystyle=4\mu+\sqrt{4}\sigma\epsilon_{4}$

The quantities $\epsilon_{i}$ ($i=1,\ldots,4$) describe the random part of the
path delays:

$\displaystyle\epsilon_{1}$
$\displaystyle=\frac{1}{2}(\xi_{12}+\xi_{24}+\xi_{46}+\xi_{67})$ (31)
$\displaystyle\epsilon_{2}$
$\displaystyle=\frac{1}{\sqrt{6}}(\xi_{12}+\xi_{24}+\xi_{43}+\xi_{35}+\xi_{56}+\xi_{67})$
$\displaystyle\epsilon_{3}$
$\displaystyle=\frac{1}{\sqrt{5}}(\xi_{12}+\xi_{24}+\xi_{45}+\xi_{56}+\xi_{67})$
$\displaystyle\epsilon_{4}$
$\displaystyle=\frac{1}{2}(\xi_{13}+\xi_{35}+\xi_{56}+\xi_{67})$

Figure 10: A pseudo-code summarising the algorithm for the calculation of the
covariance matrix of a given timing graph.

Now, considering the combinations of type
$\langle\epsilon_{i}\epsilon_{j}\rangle$ and taking into account the relations
(29), the covariance matrix (32) can be constructed. Note that it is
symmetrical since $\langle\epsilon_{i}\epsilon_{j}\rangle$ =
$\langle\epsilon_{j}\epsilon_{i}\rangle$. The procedure of the covariance
matrix estimation is summarised in the form of an algorithm, and its pseudo-
code is given in Fig. 10. The resulting covariance matrix for the graph shown
in Fig. 9 can be written as follows:

$\Sigma=\begin{pmatrix}1&\frac{3}{2\sqrt{6}}&\frac{3}{2\sqrt{5}}&\frac{1}{4}\\\
&1&\frac{4}{\sqrt{30}}&\frac{3}{2\sqrt{6}}\\\ &&1&\frac{1}{\sqrt{5}}\\\
&&&1\end{pmatrix}\,,$ (32)

where the blank spaces imply that the matrix entries are symmetric.

Figure 11: Sample timing graph with eight paths and its covariance matrix
visualised. This timing graph is used as a ‘building’ block to construct a
larger graph shown in Fig. 12. Figure 12: Cascade formed by replicating the
building block from Fig. 11 and its covariance matrix visualised. There are 64
paths in the resulting cascade.

To illustrate the algorithm and the result of its operation, we consider a
larger timing graph. Firstly, we construct a “building” block, which is a
simple timing graph containing several nodes, see Fig. 11. This simple graph
is replicated to construct a larger graph shown in Fig. 12. The application of
the algorithm shown in Fig. 10 yields the covariance matrices, also shown in
these figures. The matrices are visualised as a square field and with a colour
gradient to convey the magnitude of an entry in a matrix. We note that the
entries of a covariance matrix in are $\Sigma_{ij}\in[0,1]$, as negative
correlations between paths has no physical meaning.

##  VIII Discussion and Conclusions

In this paper, original research contributions are presented. The path-based
SSTA problem for VLSI, namely the problem of finding critical paths in a
timing graph, was mapped onto the problem of correlated extremes. For such a
problem, approximate solutions in the form of corrections to the Gumbel law
(2) for the case of weak correlations were obtained and compared with
numerical Monte Carlo simulations. In particular:

  * •

An original formulation of the problem was proposed, and the path-based SSTA
problem of a timing graph was studied from the perspective of Extreme Value
Statistics. This point of view allowed us to obtain the upper bound on the
worst case delay of VLSI. The upper bound is given by the Gumbel law.

  * •

By obtaining these corrections to the Gumbel law, the mean value of (1), which
is the total delay in a timing graph, was obtained. The total delay is
inversely proportional to the correlation coefficient $\rho$.

  * •

A general approach for studying weakly correlated Gaussian RVs, alternative to
the known techniques from Refs. [41, 42], was proposed. The starting point of
the approach was the characteristic function (11) in the case of $N$
correlated Gaussian RVs. Then, using the relaxation condition (20) for the
covariance matrix $\boldsymbol{\Sigma}$, asymptotic formulae for the CDF of
weakly correlated Gaussian RV were obtained by treating non-diagonal entries
$\varepsilon_{ij}$ of the covariance matrix as perturbations. The results (19)
and (V-A) extended the Gumbel law (2) in a form of series in powers of the
smallness parameter $\varepsilon_{ij}$.

  * •

The obtained asymptotic results were in a very good agreement with numerical
experiments. Moreover, the second order approximation (V-A) gave an
improvement over the first order one (19), as one can see from the PDFs shown
in Fig. 5. Studying the behaviour of the mean of (1) as the function of the
correlation coefficient $\rho$, it was remarkable to obtain even better
agreement with the numerical experiments for the cases of the maximum of a
relatively small number of RVs, $N=15$, $50$ and $100$. One can see from Fig.
6 that the second order approximation fits the simulations even in the regime
of strong correlations when the theory is not applicable a priori.

  * •

However, the asymptotic formulae can quickly blow for higher values of the
correlation coefficient, when the number of variables in (1) is small, as it
is shown in Fig. 7. The influences of the _topology_ of the covariation matrix
$\Sigma$ on the applicability of the corrections to Gumbel law should be
studied rigorously. For this purpose, realistic circuits should be used, which
is a scope for a separate study.

  * •

The effect of non-identical mean values and variances on the distribution of
(1) was studied numerically. It was shown that, when RVs deviate from the
identical distributed case slightly, the _average_ deviations $\delta\mu$ and
$\delta\sigma$, scaled the resulting distribution. This indicated that the
proposed theory could be applied for the non-IID case with minimal
corrections. The requirement on the deviations $\delta\mu$ and $\delta\sigma$
to be small is justified for logic circuits, as the delays in different paths
are not expected to vary greatly. The renormalisation of distributions with
respect to $\delta\mu$ and $\delta\sigma$ should be additionally studied in a
separate research.

  * •

An algorithm for estimation of the covariance matrix $\boldsymbol{\Sigma}$ was
proposed. Since it was aimed at obtaining correlations due to shared nodes,
the algorithm focused only on the random part of the delays. Assuming the
random part $\xi_{i}$ to be drawn from the standard normal distribution (30),
the combinations of the average of random components
$\langle\epsilon_{i}\epsilon_{j}\rangle$ of the path delays were calculated.
These combinations corresponded to entries of the covariance matrix.

##  IX Acknowledgement

The authors wish to thanks Paul Frain, Adrian Wrixon, Anton Belov and all the
members of the signoff team at Synopsys, Ireland, for the stimulating
discussions.

This work has emanated from research supported in part by Synopsys, Ireland,
and a research grant from Science Foundation Ireland (SFI) and is co-funded
under the European Regional Development Fund under Grant Number 13/RC/2077.

##  Appendix A On the Property of Normal Distributions

We would like to prove that the expression below is the case:

$\left.\frac{\partial\chi}{\partial\varepsilon_{ij}}\right|_{\varepsilon_{ij}=0}=\frac{1}{2}\frac{\partial^{2}\chi_{0}}{\partial\mu_{i}\partial\mu_{j}}\,.$
(33)

To prove this, compute the derivatives directly. For the characteristic
function (13), we have:

$\displaystyle\left.\frac{\partial\chi}{\partial\varepsilon_{ij}}\right|_{\varepsilon_{ij}=0}$
$\displaystyle=\left.\left[-\frac{1}{2}k_{i}k_{j}\chi_{0}(\mathbf{k})\cdot\exp\left(-\frac{1}{2}\varepsilon_{ij}k_{i}k_{j}\right)\right]\right|_{\varepsilon_{ij}=0}$
$\displaystyle=-\frac{1}{2}k_{i}k_{j}\chi_{0}(\mathbf{k}).$

The first two derivatives of the characteristic function of $N$ uncorrelated
normal RVs,
$\chi_{0}(\mathbf{k})=\exp\left(i\mu_{i}k_{i}-\frac{1}{2}\sigma_{i}^{2}k_{i}^{2}\right)$,
read

$\frac{\partial\chi_{0}}{\partial\mu_{i}}=ik_{i}\chi_{0},\quad\frac{\partial^{2}\chi_{0}}{\partial\mu_{i}\partial\mu_{j}}=-k_{i}k_{j}\chi_{0}.$

Comparing the left and right sides of the derivatives, we obtain the relation
(33). Calculating higher order derivatives, it is straightforward to show that

$\left.\frac{\partial^{2}\chi}{\partial\varepsilon_{ij}\partial\varepsilon_{kl}}\right|_{\varepsilon_{ijkl}=0}=\frac{1}{4}\frac{\partial^{2}\chi_{0}}{\partial\mu_{i}\partial\mu_{j}\partial\mu_{k}\partial\mu_{l}}.$
(34)

## References

  * [1] S. H. Gerez, _Algorithms for VLSI Design Automation_. Wiley, 1998.
  * [2] R. Ho, K. W. Mai, and M. A. Horowitz, “The future of wires,” _Proceedings of the IEEE_ , vol. 89, no. 4, pp. 490–504, Apr 2001.
  * [3] C. Visweswariah, “Death, taxes and failing chips,” in _Proc. DAC_. IEEE, Jun. 2003, pp. 343–347.
  * [4] S. Sapatnekar, _Timing_. Springer-Verlag, 2004.
  * [5] M. Orshansky, S. Nassif, and D. Boning, _Design for Manufacturability and Statistical Design: A Constructive Approach_ , ser. Series on Integrated Circuits and Systems. Springer, 2008.
  * [6] J. Bhasker and R. Chadha, _Static Timing Analysis for Nanometer Designs. A Practical Approach_. Springer, 2009.
  * [7] L. Lavagno, I. L. Markov, G. Martin, and L. K. Scheffer, Eds., _Electronic Design Automation for IC Implementation, Circuit Design, and Process Technology_. CRC Press, 2016.
  * [8] S. Wolf and R. Tauber, _Silicon Processing for the VLSI Era, Volume 1 – Process Technology_. Lattice Press, 1986\.
  * [9] M. Dietrich and J. Haase, Eds., _Process Variations and Probabilistic Integrated Circuit Design_. Springer, 2012\.
  * [10] A. B. Kahng, “New game, new goal posts: A recent history of timing closure,” in _Proc. DAC_ , June 2015, pp. 1–6.
  * [11] D. Sengupta, V. Mishra, and S. S. Sapatnekar, “Invited: Optimizing device reliability effects at the intersection of physics, circuits, and architecture,” in _Proc. DAC_ , June 2016, pp. 1–6.
  * [12] D. Saha, “Modeling process variation in scaled cmos technology,” _IEEE Design & Test of Computers_, vol. 27, no. 2, pp. 8–16, 2010.
  * [13] M. Michael, C.and Ismail, _Statistical modeling for computer-aided design of MOS VLSI circuits_ , ser. The Kluwer international series in engineering and computer science. Kluwer, 1993.
  * [14] A. Srivastava, D. Sylvester, and D. Blaauw, Eds., _Statistical Analysis and Optimization for VLSI: Timing and Power_. Springer, 2005.
  * [15] D. Blaauw, K. Chopra, A. Srivastava, and L. Scheffer, “Statistical timing analysis: From basic principles to state of the art,” _IEEE Trans. Comput.–Aided Des. Integr. Circuits Syst._ , vol. 4, no. 8, 2008.
  * [16] S. Tasiran and A. Demir, “Smart monte carlo for yield estimation,” in _Proc. TAU Workshop_ , 2006.
  * [17] R. Kanj, R. Joshi, and S. Nassif, “Mixture importance sampling and its application to the analysis of SRAM designs in the presence of rare failure events,” in _Proc. DAC_ , 2006.
  * [18] V. Veetil, D. Blaauw, and D. Sylvester, “Critically aware latin hyper- cube sampling for efficient statistical timing analysis,” in _Proc. TAU Workshop_ , 2007.
  * [19] V. Veetil, K. Chopra, D. Blaauw, and D. Sylvester, “Fast statistical static timing analysis using smart monte carlo techniques,” _IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems_ , vol. 5430, no. 6, pp. 852–865, 2011.
  * [20] R. Krishnan, D. Stamoulis, and F. Dartu, “Qcsail – quantile computation via sensitivity analysis and incremental learning,” in _Proc. TAU_ , 2018, pp. 31–36.
  * [21] H. A. David and H. N. Nagaraja, _Order Statistics_ , 3rd ed., ser. Wiley Series in Probability and Statistics. John Wiley and Sons, 2003.
  * [22] A. Lange, C. Sohrmann, R. Jancke, J. Haase, I. Lorenz, and U. Schlichtmann, “Probabilistic standard cell modeling considering non-gaussian parameters and correlations,” in _Proc. DATE_ , Mar 2014, pp. 1–4.
  * [23] A. Lange, I. Harasymiv, O. Eisenberger, F. Roger, J. Haase, and R. Minixhofer, “Towards probabilistic analog behavioral modeling,” in _Proc. ISCAS_ , May 2015, pp. 2728–2731.
  * [24] A. Lange, C. Sohrmann, R. Jancke, J. Haase, B. Cheng, A. Asenov, and U. Schlichtmann, “Multivariate modeling of variability supporting non-gaussian and correlated parameters,” _IEEE Trans. Comput.–Aided Des. Integr. Circuits Syst._ , vol. 35, no. 2, pp. 197–210, Feb 2016.
  * [25] I. Nitta, T. Shibuya, and K. Homma, “Statistical static timing analysis technology,” _Fujitsu Sci. Tech. J._ , vol. 43, no. 4, pp. 516–523, Oct. 2007.
  * [26] C. Forzan and D. Pandini, “Statistical static timing analysis: A survey,” _Integration, the VLSI Journal_ , vol. 42, no. 3, pp. 409–435, 2009, special Section on DCIS2006.
  * [27] D. Mishagli, E. Koskin, and E. Blokhina, “Path-based statistical static timing analysis for large integrated circuits in a weak correlation approximation,” in _Proc. ISCAS_ , 2019.
  * [28] C. Tracy and H. Widom, “Level-spacing distributions and the airy kernel,” _Comm. Math. Phys._ , vol. 159, pp. 151–174, 1994.
  * [29] ——, “On orthogonal and symplectic matrix ensembles,” _Comm. Math. Phys._ , vol. 177, no. 3, pp. 727–754, 1996.
  * [30] S. Majumdar and G. Schehr, “Top eigenvalue of a random matrix: Large deviations and third order phase transition,” _J. Stat. Mech. Theory Exp._ , no. 1, p. P01012, 2014.
  * [31] M. L. Mehta, _Random Matrices_ , 3rd ed., ser. Pure and Applied Mathematics. Elsevier, 2004, vol. 142\.
  * [32] F. Dyson, “Statistical theory of the energy levels of complex systems. I,” _J. Math. Phys._ , vol. 3, no. 1, pp. 140–156, 1962.
  * [33] G. Biroli, J.-P. Bouchaud, and M. Potters, “On the top eigenvalue of heavy-tailed random matrices,” _Eur. Phys. Lett._ , vol. 78, p. 10001, 2007\.
  * [34] M. Prähofer and H. Spohn, “Universal distributions for growth processes in $1+1$ dimensions and random matrices,” _Phy. Rev. Lett._ , vol. 84, no. 21, pp. 4882–4885, May 2000.
  * [35] K. Johansson, “Shape fluctuations and random matrices,” _Comm. Math. Phys._ , vol. 209, pp. 437–476, 2000.
  * [36] J. Baik and E. Rains, “Limiting distributions for a polynuclear growth model with external sources,” _J. Stat. Phys._ , vol. 100, pp. 523–541, 2000.
  * [37] G. Lemarié, A. Kamlapure, D. Bucheli, L. Benfatto, J. Lorenzana, G. Seibold, S. C. Ganguli, P. Raychaudhuri, and C. Castellani, “Universal scaling of the order-parameter distribution in strongly disordered superconductors,” _Phys. Rev. B_ , vol. 87, no. 18, p. 184509, 2013.
  * [38] M. Evans and T. Hanney, “Nonequilibrium statistical mechanics of the zero-range process and related models,” _J. Phys. A: Math. Gen._ , vol. 38, no. 19, pp. R195–R240, 2005.
  * [39] J.-Y. Fortin and M. Clusel, “Applications of extreme value statistics in physics,” _J. Phys. A: Math. Theor._ , vol. 48, p. 183001, 2015.
  * [40] S. N. Majumdar, A. Pal, and G. Schehr, “Extreme value statistics of correlated random variables: A pedagogical review,” _Physics Reports_ , vol. 840, pp. 1–32, Jan. 2020.
  * [41] G. Györgyi, N. R. Moloney, K. Ozogány, and Z. Rácz, “Finite-size scaling in extreme statistics,” _Phyr. Rev. Lett._ , vol. 100, no. 21, p. 210601, 2008.
  * [42] G. Györgyi, N. R. Moloney, K. Ozogány, Z. Rácz, and M. Droz, “Renormalization-group theory for finite-size scaling in extreme statistics,” _Phys. Rev. E_ , vol. 81, no. 4, p. 041135, 2010.