††thanks: 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.

# Gate–Level Statistical Timing Analysis:
Exact Solutions, Approximations and Algorithms

Dmytro Mishagli<EMAIL_ADDRESS>Eugene Koskin  Elena Blokhina
<EMAIL_ADDRESS>School of Electrical and Electronic Engineering,
University College Dublin, Belfield, Dublin 4,
Dublin, Ireland

###### Abstract

In this paper, the Statistical Static Timing Analysis (SSTA) is considered
within the block–based approach. The statistical model of the logic gate delay
propagation is systematically studied and the exact analytical solution is
obtained, which is strongly non-Gaussian. The procedure of handling such (non-
Gaussian) distributions is described and the corresponding algorithm for the
critical path delay is outlined. Finally, the proposed approach is tested and
compared with Monte Carlo simulations.

Statistical Static Timing Analysis; non-Gaussian; delay; Integrated Circuit

###### Contents

  1. I Introduction
  2. II Mathematical Preliminaries
  3. III Background and State of the Art
    1. III.0.1 Phenomenological delay models
    2. III.0.2 Parameterised delay models
    3. III.0.3 Microscopic delay models
  4. IV Statement of the Problem
  5. V Gate Delay Distribution
    1. V.1 Maximum of Two Correlated Gaussians
    2. V.2 Derivation of the Exact Formula
    3. V.3 Limit Cases of the Correlation Coefficient
  6. VI Moments of the Gate Delay Distribution
    1. VI.1 Mean
    2. VI.2 Standard Deviation
    3. VI.3 Higher Moments
  7. VII Delay Propagation Algorithm
    1. VII.1 Brief Theoretical Summary
    2. VII.2 Algorithm
    3. VII.3 Optimisation
  8. VIII Discussion and Conclusions

##  I Introduction

Timing verification has never been more important for the design of digital
integrated circuits (ICs) than now. Indeed, fabricated circuits must safely
operate under working conditions and should satisfy the designed timing
constraints. With the decrease of feature size, the gap between designs and
fabricated circuits widens significantly. This is not only due to the
increasing complexity of designs, but also because of variations of the
parameters. At the moment, the technology has already scaled down below 5 nm
getting close to physical limits. At such scales, the occurring variations of
feature parameters (e.g. channel length) are of order of their nominal values.
The traditional methods of the timing verification, static timing analysis
(STA) and Monte Carlo (MC) simulations, no longer give satisfactory results.
STA has become too pessimistic as it relies on utilising so-called corner
values of the parameters, i.e. takes the worst case scenario ($3\sigma$
deviation from nominal values) even if the probability of it to happen is
negligible in practice. At the same time, MC, while remains the most reliable
tool, can take enormous time to finish a single design. All this moves
research within statistical static timing analysis (SSTA) to front.

The main advantage of SSTA (and its curse at the same time) is that all
design’s parameters are considered as random variables (RVs) from the very
beginning. It is especially important now as the impact of truly random
variations (such as dopant concentration or oxide thickness) has become
significant. Thus, the use of statistical approaches is natural. The main goal
and task of SSTA is to determine a final distribution, probability density
function (PDF) or cumulative density function (CDF), of a design slack, a so-
called critical path delay.

From the analysis point of view, an IC is a so-called _timing graph_ [1].
Hence, the problem of timing analysis (either STA or SSTA) is equivalent to a
graph optimisation problem, where one is required to find a shortest (longest)
path in a graph [2]. There are two general approaches to the problem:

  * (i)

block-based methods, which consider a delay propagation iteratively, level by
level, from source to sink;

  * (ii)

path-based methods, within which a timing graph is considered as a set of
paths.

For the reasons that will become clear further in this paper, the most of
researchers’ activity was carried out within the block-based approach.

Figure 1: Combinational Logic Circuit modelling in Static Timing Analysis. A
simple logic circuit on the left, and its timing graph in the middle and on
the right. The timing graph can be traversed either within block-based or
path-based approach.

The statistical domain brings the problem to another level of complexity,
graph optimisation under uncertainties. This rises new challenges, such as
computing a maximum of correlated non-Gaussian RVs. These challenges as well
as various attempts to mitigate them have been described in two excellent
reviews by Blaauw, Chopra et al [3] and Forzan and Pandini [4], devoted to the
first decade of SSTA (and mainly to block-based approaches). One of such
challenges, which is of current interest, is the non-Gaussian nature of
delays. Logic gate operation time is determined by the _latest_ of its arrival
signals. In other words, nodal delay calculation requires computation
$\max(\ldots)$ of all arrival signals’ times. Even assuming the signals to
have Gaussian delay distributions, one faces a problem that the result of
maximum operation is a non-Gaussian PDF. Traditionally, at this stage
simplifications and$/$or approximations are introduced that ineluctably bring
errors to the delay computation. Moreover, it is widely believed in the
community that exact solution does not exist (see e.g. [5]).

This paper concludes our research started in [6, 7]. For the Gaussian arrival
delays, an exact expression of the distribution of a logic gate delay
neglecting correlations was presented without derivation in [6]. In principle,
any non-Gaussian PDF can be represented as a mixture of Gaussian kernel
functions, which allows one to apply the expression and propagate the delay
with almost no loss in information. Such a representation requires solving of
a corresponding optimisation problem. This was then investigated in [7], where
it was shown that a concept of a _Gaussian comb_ can give error in mean
compare to MC less than $0.01\%$. In this paper, we present derivation of the
logic gate delay distribution in a general case, taking correlations between
nodes into account. We also investigate first four moments of the
distribution: mean, standard deviation, skewness and kurtosis. Having full
mathematical description of a system, we then discuss graph traversal
algorithm that includes solution to an optimisation problem. Finally, we
verify our approach by considering several application scenarios.

The paper is organised as follows. In Section III we give details on the
background and highlight the most common approaches and latest results in the
field. The general statement of the problem is given in Section IV. The core
of the paper is in the next sections. Hence, we give the derivation of the
gate delay distribution in Section V. The algorithm for decomposition of non-
Gaussian distributions and timing graph traversal is presented in Section VII.
Section VIII gives discussion of the obtained results.

##  II Mathematical Preliminaries

In this work we will use the following explicit notation for the Gaussian PDF:

$\phi(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}},$
(1)

where the function $\varphi(x)$ is referred to as the _Gaussian kernel_. The
CDF of the _standard_ Gaussian distribution will be denoted as follows:

$\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],$
(2)

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}.$ (3)

For a compact notation, we will use

$F_{i}(x)\stackrel{{\scriptstyle\text{def}}}{{=}}\Phi\left(\frac{x-\mu_{i}}{\sigma_{i}}\right)=\frac{1}{2}\left[1+\operatorname{erf}\left(\frac{x-\mu_{i}}{\sqrt{2}\sigma_{i}}\right)\right].$
(4)

Let the function $g(z)$ be the PDF of the maximum
$\zeta_{2}=\max(X_{1},X_{2})$ for two independent RVs $X_{1}$ and $X_{2}$
distributed normally. This PDF can be easily found. In this case, we note that
the PDF and CDF can be factorised: $f(x_{1},x_{2})=f(x_{1})f(x_{2})$ and
$F(x_{1},x_{2})=F(x_{1})F(x_{2})$. In this case, we have:

$g(z)=f_{1}(z)F_{2}(z)+f_{2}(z)F_{1}(z),$ (5)

where $f_{i}(x)$ and $F_{i}(x)$ are a Gaussian PDF and CDF of the
$i^{\text{th}}$ distribution ($i=1,2$) respectively.

The PDF for the maximum of two _correlated_ Gaussian RVs has the form:

$\displaystyle g(z,\rho)$
$\displaystyle=f_{1}(z)\Phi\left[\frac{1}{\sqrt{1-\rho^{2}}}\left(\frac{z-\mu_{2}}{\sigma_{2}}-\rho\frac{z-\mu_{1}}{\sigma_{1}}\right)\right]$
(6)
$\displaystyle+f_{2}(z)\Phi\left[\frac{1}{\sqrt{1-\rho^{2}}}\left(\frac{z-\mu_{1}}{\sigma_{1}}-\rho\frac{z-\mu_{2}}{\sigma_{2}}\right)\right].$

##  III Background and State of the Art

In this section, we briefly outline the principles of SSTA and overview
traditional techniques as well as recent papers.

In STA, whether it is deterministic or statistical, a combinational logic
circuit is modelled with a timing graph, as it is shown in Fig. 1. The graph
is a Direct Acyclic Graph, i.e. it is given by a set $V$ of vertices and a set
$E$ of edges, $G(V,E)$. Vertices correspond to nodes (logic gates) while edges
describe interconnect delays. As we noted in Introduction, there are two
methods to traverse the timing graph, path- and block-based. Within the path-
based method, the total circuit delay equals to the maximal path delay,
namely, for a circuit with $n$ paths, we have

$D=\max(\tau_{1},\tau_{2},\ldots,\tau_{n}),$ (7)

where $\tau_{i}$ is an accumulated delay in a path $i$. In SSTA, all
$\tau_{i}$ are RVs and, in general case, are correlated and non-Gaussian. The
distribution of the circuit delay $D$ in such a case is a non-trivial problem
and remains unsolved. We will address this issue in a separate study.

Within the block-based method, the graph is traversed in a levelised manner,
i.e. the delay is propagated from source to sink of a graph and computed for
each block (see Fig. 1). One cannot propagate further in a graph until all
delays are computed within a block. The maximum transition time is determined
by the latest arrival signal, thus, for a gate with two inputs, we have

$\ldots+\max(\tau_{1},\tau_{2})+\tau_{\text{g}}+\tau_{\text{int}}+\ldots,$ (8)

where $\tau_{1}$ and $\tau_{2}$ are the arrival times, $\tau_{\text{g}}$ is a
gate’s characteristic operation time, $\tau_{\text{int}}$ is an interconnect
delay. All other delays like a wire delay $\tau_{\text{w}}$ can be added to
(8) and the sum can be continued. We will represent an output of a logic gate
as follows:

$D_{\text{gate}}=\max(\tau_{1},\tau_{2})+\tau_{0},$ (9)

where $\tau_{0}$ is a delay associated with interconnects, wires, etc. Thus,
the block-based approach requires computation of (9) at each node of a graph.
In SSTA $\tau_{i}$ are RVs, and this expression is a source of all challenges
of the approach.

Most of the works on SSTA are performed with the block–based approaches as
attacking the SSTA problem (LABEL:chap3:eq:ssta_problem) is computationally
hard because the total number of paths in IC grows exponentially. The research
is focused on developing models of delays with further incorporating the
models into algorithms for SSTA. One can distinguish three categories of the
delay models:

  1. (i)

Phenomenological, or non-parameterised, models. The approaches based on these
models assume distributions of delays to be given (e.g. from pre-
characterisation of circuits elements using SPICE) and propose a framework to
propagate the distributions through the timing graph $G$.

  2. (ii)

Parameterised models. These approaches model delays as functions of PVT
variations and work directly with RVs. Information regarding the variations
can be obtain from foundries or via MC simulations. As their final step, the
methods need to relate the models to particular distributions. Techniques like
AWE [8, 9] can be used for this.

  3. (iii)

Microscopic models aim to establish relations between the PVT variations and
delays in transistors and, hence, logic gates. These models are developed for
better understanding of the problem and for more accurate parameterised models

Below we review these models in detail.

####  III.0.1 Phenomenological delay models

A pioneering work on SSTA was done by Berkelaar [10], where it was proposed to
consider distributions of the arrival times’ delays of a logic gate, i.e.
$\tau_{1}$ and $\tau_{2}$ in (LABEL:chap3:eq:gate_delay_2), as independent
Gaussian RVs. Thus, the distribution of $\max(\tau_{1},\tau_{2})$ is given by
(5). It was proposed then to approximate the mean value of a gate delay simply
by adding the mean value of $\max(\tau_{1},\tau_{2})$ to the mean value of
$\tau_{0}$, thus, treating the distribution (5) as if it was a Gaussian one.
The standard deviation of a gate delay is calculated as
$\sigma^{2}_{\text{gate}}=\sigma^{2}_{\max}+\sigma^{2}_{0}$. After these
calculations performed, the gate delay was considered as a Gaussian RV, i.e.
its distribution given by (1), and the whole process was repeated at the next
gate. In principle, this repeats the famous Clark’s algorithm [11], but for a
case of zero correlations. The exact expression for (5), as well as for the
mean value and standard deviation of an RV given by this distribution, were
presented later, in a joint work with Jacobs in [12], where this delay model
was applied to a gate sizing problem.

In these works, [10] and [12], correlations between delays were not taken into
account. Tsukiyama et al. addressed this issue in [13, 14]. It was proposed to
use the known expressions for the PDF of two correlated Gaussian RVs, namely,
equation (6). Other steps were the same as in works by Berkelaar and Jacobs.
Thus, for each gate (a node in a timing graph), the exact distribution (6) was
approximated with a Gaussian distribution by matching two first moments and
taking into account the correlation coefficient, as per Clark [11]. Such an
algorithm has a complexity of $\mathcal{O}(|V|\cdot|E|)$, as it was reported
in [14]. Hence, the approach by Berkelaar [10] has $\mathcal{O}(|V|+|E|)$ time
complexity, as the delays are independent.

Azuma et al. [15] have proposed to approximate the PDF of the maximum of two
Gaussian distributions via GMM. The authors gave closed formulas for the
parameters of the GMM. It should be noted that the proposed approach required
that the distributions one takes the maximum of are either Gaussians or given
in the GMM form. However, it needs to be explained how an arbitrary non-
Gaussian distribution should be decomposed into the mixture. Therefore,
assuming the initial delays are distributed normally, the overall complexity
with respect to the total number of atomic operations, $N$, is
$\mathcal{O}(m^{2}N)$, where $m$ is the number of components in the GMM.

Chen et al. [16] proposed to approximate non-Gaussian distributions by an
inverse Gaussian distribution, which has a shape parameter responsible for
skewness. Such an approximation is meant to be done as a preprocessing step,
while the actual distributions can be retrieved from prior simulations (this
can be done only ones). Then, the traversal of these non-Gaussian
distributions can be done (hence, convolution) via summation of the parameters
of the inverse Gaussian PDFs. The proposed approach works well when the $\max$
is not required, e.g. for a chain of inverters of for a cascade of AND (NAND)
gates, when both inputs are equal (thus, $\max$ is trivial). This makes the
approach not suitable for real VLSI designs.

A very recent work by Jin et al. [17] uses a so-called
log–extended–skew–normal distribution to represent gate delays. This
distribution has four parameters to match not only the mean and standard
deviation, but also the skewness and kurtosis. Though it has remarkable
agreement with MC simulations, it has the same limitations as the work by Chen
et al. [16], i.e. it is unclear how the $\max$ operation should be performed.

####  III.0.2 Parameterised delay models

Seminal works within the block–based approach were made by Chang & Sapatnekar
[18, 19], Chang et al. [20], and Visweswariah et al. [21, 22]. These works
continued the approach of Berkelaar to use the moment matching Clark’s method
[11] and have shaped block–based SSTA. The works by Visweswariah et al. [21,
22] resulted in an IBM’s tool EinsStat for SSTA. Let us consider these papers
in detail.

It was proposed to consider a gate delay as a function of process variations.
Thus, for $n$ sources of variations $X_{i}$, the gate delay can be written in
a general form as

$d=d(X_{1},X_{2},\ldots,X_{n}),$ (10a) where the exact form of the r.h.s.
depends on a delay model used. Thus, Chang & Sapatnekar [18, 19] proposed a
direct extension of [10, 12, 13, 14] by transforming a set of correlated RVs
into a set of i.i.d. that are distributed normally, i.e.
$\sim\mathcal{N}(0,1)$, by PCA. In this case, (10a) is a linear combination of
variation parameters. For such linear combinations, the Clark’s algorithm [11]
is then used.

A notable result of [18, 19] is a method to extract correlations. Chang &
Sapatnekar proposed a grid model for spatial correlations, where a die is to
be split into $n$ grids, and the correlation between circuit elements that are
in the grids decays with the distance between the grids. The overall
complexity of the algorithm is $\mathcal{O}(pn(|E|+|V|))$, where $p$ is a
number of variation parameters considered and $n$ is a number of grids. It
also should be noted that a quad–tree method for modelling the spatial
correlations on a die was proposed by Agarwal et al. [23]. As the name
suggests, the die is divided into levels, and for each level $i$, the die is
partitioned into $2^{(i)}$-by-$2^{(i)}$ squares, where the square for the
first level has an area of the whole die, and the last level is split into 16
squares. An independent RV is then associated with each square. These grid and
quad–tree methods are often employed in different parameterised block–based
approaches. Further information can be found in a review by Forzan & Pandini
[4].

At the same time, Visweswariah et al. [21, 22] proposed a linear _canonical
model_ of a delay. Thus, the general delay (10a) was written as follows:

$d=a_{0}+\sum\limits_{i=1}^{n}a_{i}\Delta X_{i}+a_{n+1}\Delta R,$ (10b)

where $a_{0}$ is the mean or nominal value, $\Delta X_{i}$ represent the
variation of $n$ global sources of variation $X_{i}$ from their nominal
values, $a_{i}$ are sensitivities to each of the RVs, and $\Delta R$ is the
variation of an independent RV $R$. The sensitivities can be obtained, for
example, from circuit simulations when a chip is characterised. The Clark’s
method [11] was written in terms of probabilities, using a concept of
_tightness probability_ , $T_{A}=P(A>B)$, a probability that arrival time $A$
is greater than $B$. Thus, closed–form expressions for a linear approximation
of two canonical forms (10b) was obtained. At this point, these works are
similar to [18, 19], with the only difference that the variance of the
approximation of two linear forms is computed differently. This approach has
linear complexity with a circuit size, $\mathcal{O}(n\cdot N)$, where $N$ is
the total number of $\max$ and $\mathrm{sum}$ operations.

The linear canonical model (10b) used in the discussed above approaches, is
one of the main sources of inaccuracies, which was quickly understood. This is
because the linearity of (10b) implies Gaussian distribution for delays and
the $\max$ operation. Therefore, many attempted to address this problem,
remaining in a paradigm of the canonical forms (10a). Thus, Chang et al. [20]
proposed the _generalised canonical form_

$d=a_{0}+\sum\limits_{i=1}^{n}a_{i}\Delta X_{i}+a_{n+1}\Delta R+f_{d}(\Delta
X_{N}),$ (10c)

where $f_{d}(\Delta X_{N})$ represents dependencies of $d$ on nonlinear
(hence, non-Gaussian) parameters. Chang et al. proposed to use numerically
computed tables to describe $f_{d}(\Delta X_{N})$.

Alternative to adding non-linear terms to the canonical form approach was
proposed by Singh and Sapatnekar [24, 25]. The delay was still represented
using the linear canonical model (10a) where non-Gaussian RVs were given by
another term. The non-Gaussian term is then addressed by _preprocessing_ the
random vector of correlated random variables via [26, 27, 28]. The correlated
Gaussian part of (10a) is handled with PCA as in [18, 19]. After the
preprocessing step, which is supposed to be made only once, hence, it does not
contribute to the total complexity of the algorithm, the moments of the
transformed RVs are calculated. The PDFs and CDFs are restored from the
canonical forms using numerical techniques proposed in [29]; the techniques
uses Padé Approximants and requires computation of moments of the canonical
forms.

In contrast to [20], Zhan et al. [30] and Zhang et al. [31] introduced
quadratic terms into (10b) instead of a function $f_{D}(\Delta X_{N})$.
However, such an approximation of the $\max$ still introduces noticeable
errors. For $N$ being a circuit size, the worst–case timing complexity of both
approaches is $\mathcal{O}(n^{2}N)$, where $n$ is a number of terms in the
quadratic canonical form. The authors indicated that the complexity can be
reduced to $\mathcal{O}(n\cdot N)$ if the most of the cross terms ignored. At
the same time, [30] uses expensive numerical integration for computation of
the $max$ of the canonical forms, and [31] requires additional
simplifications. Therefore, these approaches do not give advantage compare to
the linear canonical model.

Reference |  Complexity | Features |  Comments
---|---|---|---
[10, 12] |  $\mathcal{O}(|V|+|E|)\phantom{\displaystyle\frac{1}{2}}$ | G., $\rho=0$ |  Eq. (5) is approximated with a Gaussian by two moments matching
[13, 14] |  $\mathcal{O}(|V|\cdot|E|)$ | G., $\rho\neq 0$ |  Correlations are taken into account by utilising Eq. (6); the result is approximated with a Gaussian by matching two moments.
[18, 19] |  $\mathcal{O}(pn\cdot N)$ | G., $\rho\neq 0$ |  Correlations are treated with PCA; $\max$ is approximated with a Gaussian by Clark’s method [11];
$p$ is a number of variation parameters, $n$ is a number of grids.
[20, 21, 22] |  $\mathcal{O}(n\cdot N)$ | G., $\rho\neq 0$ |  Linear canonical form is used and, hence, $\max$ is approximated with a Gaussian using tightness probability; correlations are due to global sources of variations, $n$.
[24, 25] |  $\mathcal{O}(n\cdot|E|)$ | non-G., $\rho\neq 0$ |  Linear canonical model is used, but non-Gaussian RVs are treated with ICA. Correlated RVs are transformed with PCA as in [18, 19]; pdfs are restored from the canonical forms using Padé-Approximants-based methods [29].
[30] |  $\mathcal{O}(n^{2}N)$
$\approx\mathcal{O}(n\cdot N)$
| non-G., $\rho\neq 0$ |  Quadratic canonical form is used; $\max$ is computed using numerical techniques; Correlations are taken effectively, as in [20, 21, 22];
$n$ is the same as above.
[31] |  $\mathcal{O}(n^{2}N)$
$\approx\mathcal{O}(n\cdot N)$
| non-G., $\rho\neq 0$ |  Quadratic canonical form is used; Correlations and $\max$ are treated in the same manner as in [18, 19, 20, 21, 22];
$n$ is the same as above.
[32] |  $\mathcal{O}(n^{2}N)$ | non-G., $\rho\neq 0$ |  Delays are represented via Taylor series; $\max$ is represented via polynomials by linear regression; further matching of the three moments gives a new canonical form;
$n$ is the highest order of the polynomials.
Table 1: Complexity of block–based SSTA approaches from the literature. Here $|E|$ is a number of edges, and $|V|$ is a number of vertices in a timing graph; $N$ is a total number of $\max$ and sum operations. Further details are in the text. Reference |  Complexity | Features |  Comments
---|---|---|---
[33, 34] |  $\mathcal{O}(nK^{3}N)\phantom{\displaystyle\frac{1}{2}}$
$\approx\mathcal{O}(nK^{2}\cdot N)$
| non-G., $\rho\neq 0$ |  PDF of $\max$ of two quadratic canonical forms is approximated by _Fourier series_ ; further matching of the three moments gives a new canonical form;
$K$ is the highest order of Fourier series, $n$ is a number of variation
sources.
[35, 36] |  $\mathcal{O}(n^{3}N)$
$\approx\mathcal{O}(n\cdot N)$
| non-G., $\rho\neq 0$ |  Extension of [33, 34]; here PDF of $\max$ of two quadratic canonical forms is approximated by _least–squares–error minimisation_ ;
$n$ is as above.
[37, 38] |  $\mathcal{O}(n^{2}N)$ | non-G., $\rho\neq 0$ |  PDF of $\max$ and resulting delays are modelled with skew–normal distribution; delays are propagated via moment matching;
$n$ is as above.
[39] |  $\mathcal{O}(m^{2}N)$ | non-G., $\rho=0$ |  Delays are represented via moments, and atomic operations are applied to them; PDF is restored via AWE method;
$m$ is the number of moments used.
[16] |  NA | non-G., $\rho=0$ |  Non-G. is approximated by an inverse Gaussian distribution, and the convolution is done via summation of the parameters of such distributions;
$\max$ _is not addressed_.
[15] |  $\mathcal{O}(m^{2}N)$ | G., $\rho\neq 0$ |  $\max$ of two Gaussians is approximated with a GMM with further model reduction;
$m$ is the number of components in GMM.
[5] |  $\mathcal{O}(nm\cdot N)$ | G., $\rho\neq 0$ |  Linear canonical delay model is used; $\max$ is performed as pairwise comparison of polynomials;
$n$, $m$ are a total numbers of polynomials in a fanin set for a node and in
its delay set respectively.
Table 2: Complexity of block–based SSTA approaches from the literature
(continuation of Table 1). Here $N$ is a total number of $\max$ and sum
operations. With “NA” marked those approaches where the total complexity is
unknown since the authors did not propose the full SSTA algorithm. Further
details are in the text.

A different approach was proposed by Khandelwal and Srivastava [32], where the
authors did not rely on the canonical form (10) explicitly. Instead, gate
delays and arrival times were presented as polynomials using a Taylor series
expansion. With the second order degree, this approach is similar to (10c)
though. The maximum of two polynomials was approximated with a polynomial as
well by linear regression, instead of the moment matching technique. However,
the regression requires a MC simulation in the inner loop of the proposed
procedure. The spatial correlations are encountered effectively as
coefficients of the polynomials with the grid–based correlation model
employed. Thus, the proposed method scales linearly with the total number of
$\max$ and $\mathrm{sum}$ operations, and it is quadratic with respect to the
degree of polynomials used, $\mathcal{O}(n^{2}N)$. But the hidden MC
simulation makes this approach unlikely to be used for real VLSI designs.

The above discussed approaches either were inaccurate (since rely on a linear
model of a delay) or did not scale well due to usage of expensive numerical
methods. Therefore, Cheng et al. attempted to mitigate these disadvantages in
a series of works, [33, 34] and [35, 36]. The first approach [33, 34] proposes
another quadratic canonical model of delay:

$D=d_{0}+\sum_{i}^{n}(a_{i}X_{i}+b_{i}X_{i}^{2})+a_{\text{r}}X_{\text{r}}+b_{\text{r}}X_{\text{r}}^{2},$
(10d)

where $d_{0}$ is nominal delay (mean value), $X_{i}$ is an RV that represents
global sources of variation and $X_{\text{r}}$ is a purely independent random
variation. The RVs $X_{i}$ allowed to follow arbitrary distribution, and the
correlations are proposed to be addressed via ICA. The bottleneck of the
method is the PDF of the $\max$ operation applied to the canonical forms
(10d), which is computed using Fourier series. Thus, the overall complexity of
this method with respect to a total number of atomic operations, $N$, is
$\mathcal{O}(nK^{3}N)$, where $n$ is a total number of variation sources and
$K$ is the highest order of Fourier series. Then, the complexity was decreased
to $\mathcal{O}(nK^{2}N)$ in [36].

The second study of Cheng et al. [35, 36] goes further towards reducing the
complexity of [33]. It was proposed a different and more efficient way of
constructing PDFs from the maximum of two canonical forms. Similar to [40],
the non-Gaussian RVs were modelled with a quadratic polynomials. The $\max$
was approximated with another quadratic polynomial, $h(A,B,\Theta)$, where
$A$, $B$ are non-Gaussian RVs written in a quadratic canonical form and
$\Theta$ is a set of three fitting parameters. Then, these parameters were
found from an optimisation problem of minimising a quadratic error between
$h(A,B,\Theta)$ and $\max(A-B,0)$, while matching the mean value of the $\max$
operation. There was proposed a closed–from solution to the optimisation
problem, thus, this gave a significant speed up compare to [33, 34]. For the
total number of $n$ variational sources, the overall complexity of this method
is $\mathcal{O}(n^{3}N)$, however, it can be reduced to $\mathcal{O}(n\cdot
N)$ if the cross–terms in the canonical forms are ignored similar to [30, 31].

A different approach was initiated by Chopra et al. [41], who proposed to take
into account skewness of the $\max$ operation by modelling it with a
skew–normal distribution [42, 43], which PDF is
$f(x;\lambda)=2\phi(x)\Phi(\lambda x)$, where $\lambda$ is a shape parameter
that determines the skewness. The authors approximated the $\max$ of two
linear canonical forms (10b) with the skew–normal distribution by
moment–matching techniques, similar to Clark’s [11], and derived corresponding
formulas to take into account the third moment.

This direction was continued in the works by Vijaykumar & Vasudevan [37] and
Ramprasath et al. [38], who extended this idea beyond the linear canonical
form. Thus, the authors of [38] introduced the so-called skew-canonical
representations of the quadratic canonical form (10d) and its linear
relaxation. The moments of the skew–normal approximation of the $\max$ are
expressed via Owen’s $T$-function [44]. The correlations were represented
using the quad–tree model [45]. The authors reported that the accuracy of
their approach is comparable to that of Cheng’s et al. [35, 36] but the run
time on average was two orders of magnitude lower than for quadratic models.
This is due to the fact that this approach does not require the restoration of
a PDF at each step neither by Fourier series nor by least–squares fitting. One
can conclude that the complexity of the approach is comparable to quadratic
models with moment matching approaches, $\mathcal{O}(n^{2}N)$.

Shebaita et al. [39] proposed to represent delays via moments in terms of
distributions’ poles and residues. The authors gave general formulas for the
moments given an arbitrary distribution. Obtaining the poles and residues of
the distributions is the required pre-prepossessing step for the method. The
resulting distributions, PDFs and CDFs, are proposed to be restored using the
AWE method [8, 9]. However, the problem of correlations is not addressed in
the paper, which limits the applicability of the approach. Given the number of
moments used to represent delays, $m$, the overall timing complexity with
respect to the total number of atomic operations, $N$, is
$\mathcal{O}(m^{2}\cdot N)$.

A recent work by Lewis & Schmit [5] aims to decrease the number of RVs used in
parameterised block–based approaches. Using the linear canonical model (10a)
as a starting point, the delays are represented in terms of polynomial basis
functions, which are then propagated through a timing graph. The information
about spatial correlations is embedded into the polynomials from grid–based
models of the correlations. The $\max$ is performed as pairwise comparison of
the polynomials, and the summation results in merging the polynomials in a new
one. However, the proposed methodology is limited only to Gaussian RVs, and
the authors did not discuss possible extension to non-Gaussian distributions.
The overall complexity with respect the total number of atomic operations,
$N$, is $\mathcal{O}(nm\cdot N)$, where $n$ is a total number of polynomials
in the fanin set for a node, and $m$ is a total number of polynomials in the
delay set for a node.

Another recent work by Jin et al. [46] aims to increase the accuracy of
block–based approaches by treating the $\max$ operation using multivariate
adaptive regression splines (MARS). The idea is to use SPICE simulations to
obtain training data to build the MARS models. The result of the $\max$
operation is mapped onto a skew–normal distribution. The complexity of the
algorithm with respect to the number of the atomic operations, $N$, is
$\mathcal{O}(MN)$, where $M$ is the number of samples in the training data.
Overall, while this approach has high accuracy, it requires exhaustive
numerical computations at each step.

To summarise, we have seen that all the proposed methods scale linearly with
the total number of the atomic ($\max$ and $\mathrm{sum}$) operations in a
circuit, $N$, but with a scale factor, i.e. $\mathcal{O}(k\cdot N)$ (see
Tables 1 and 2). This scale factor, $k$, varies with approaches and can imply
exhaustive numerical computations at each gate. Linearity with $N$ should not
be surprising as all the methods use BFS algorithm, which complexity is
$\mathcal{O}(|V|+|E|)\approx\mathcal{O}(N)$.

####  III.0.3 Microscopic delay models

There is a brunch of research dedicated to studying the sources of variations,
which is of particular importance for effective delay models. The main focus
is on developing the current–based models of delays of the logic cells [47].
Since these works aim to determine and$/$or model the fundamental stochastic
behaviour of the process variations, we group them as _microscopic_ models,
even though they do not consider explicitly the microscopic properties of
matter such as distributions of particles. Developing a microscopic model was
not aimed to be a part of the Thesis, thus, we will only briefly discuss some
of the works.

Thus, Fatemi et al. [48] modelled the cell output voltage waveform as a
Markovian stochastic process considering the variations of physical parameters
(such as channel width, length, and temperature) Gaussian. This allowed the
authors to obtain a distribution of the cell delay, which was in better
agreement with the MC simulations done in SPICE than the most recent
current–based model at that time by Keller et al. [49].

In parallel, Y. Cao & L.T. Clark [50] presented a physically based, analytical
delay variability model for both regions, at saturation and at sub-threshold.
Starting from the well-known $\alpha$-power law for a gate delay in a MOSFET
[51], the authors took short–channel effects into account based on the
previous work by Orshansky et al. [52]. That work allowed the authors to (i)
model the dependence of $V_{\text{th}}$ on channel length, (ii) utilise a
formula for the saturation of carrier velocity, and (iii) formulate a unified
formula for the drive current.

Shinkai et al. [53, 54] considered a wide range of PVT variations and their
influence on the gate delay. With the focus on fluctuations of the output
current, $I_{\text{d}}$ the authors proposed a methodology to relate PVT
variations to the gate delay and output slew given input slew and output load.
It should be noted that the methodology is agnostic with respect to the gate
delay model, thus, any of the models discussed in the previous subsection can
be applied.

Sinha et al. [55] studied the possibility of practical usage of the
current–source models in industrial timing analysis environments [56]. The
authors have presented a technique for efficient storing the waveforms in
memory and validated it for 14 nm technology designs.

P. Cao et al.[57] related analytical results from [58, 59] for the current
variations at sub- and near-threshold regions to the parameters of the
log–skew–normal (LSN) distribution by moment matching. In other words, they
gave microscopic reasoning for gate delays modelled with LSN distributions.

##  IV Statement of the Problem in Block – Based SSTA

Consider a simple combinational logic circuit and its timing graph, as shown
in Figure 2. As was noted above, the edges of this graph represent logic
gates, and the vertices represent the inputs and outputs of the logic gates.
Firstly, we note that logic gates have an ‘internal structure’ since they are
made of transistors. This results in a finite characteristic time needed for
gates to operate. This is one of the sources of delays in a circuit. Secondly,
due to operational delays existing in every gate (and also in every
interconnect), signals, propagating in a logic circuit, arrive at the next
logic gate also with delays. Most logic gates have at least two inputs. They
need both inputs to ‘arrive’ (i.e. to change from Low to High or vice versa)
before a logic operation can be performed on them and the output can be
generated. For this reason, the maximum of two (or more) input delays is
another contribution to circuit delay. We will discuss both contributions and
their mathematical form below.

Thus, the main problem of block-based timing analysis for logic circuits can
be formulated as the mathematical problem of calculating the $\max$ function
of arrival times and accommodation the operation time of logic gates and
interconnect in the overall delay calculation. In other words, this is a
problem for graph optimisation. Graphs are typically described by
corresponding adjacency matrices. Firstly, the BFS algorithm (see Figure
LABEL:chap2:fig:pseudo-code_breadth) is used to understand in what order a
timing graph is traversed. As soon as the sequence of operations required to
find the delay between two nodes is found, one must use gate delays and
interconnect delays involved in this sequence to calculate the actual delay.

Both contributions, the gate operation time and the $\max$ function, can have
a significant impact on the overall circuit delay. It is difficult to say
which contribution is more important, and they both should be taken into
account. In principle, it is quite straightforward to do in the context of
deterministic timing analysis, but it is not the case when uncertainty arises.
When it is necessary to include variations of parameters, SSTA is required. In
such a case, the arrival and gate operation times are described by random
variables given by some distributions.

Figure 2: An example combinational circuit with a symbolic notation of logic
gates and interconnects. The timing graph is split in a levelised manner
(blocks), which indicates the block–based analysis.

Looking at the mathematical statement of the problem, we recall that at the
individual gate level, delay propagation is described by two _atomic_
operations: computing the maximum ($\max$) of two delays entering a gate and
the summation ($\mathrm{sum}$) of the latter with the delay of the gate due to
its operation time, the interconnect delay, the wire delay, etc. From the
statistical point of view, when these operations are applied to RVs, the delay
of a gate with two inputs reads

$\eta=\max(X_{1},X_{2})+X_{0},$ (11)

where $X_{1}$ and $X_{2}$ are the RVs that describe the arrival times of input
signals, and $X_{0}$ is the RV that gives the gate operation time.

In terms of distributions, combinations like (11) are equivalent to the
_convolution_ of two functions: the PDF $f_{\text{max}}(x)$ of the maximum
$\max(X_{1},X_{2})$ of the arrival signals’ delays with the PDF $f_{0}(x)$
that gives the distribution of the gate operation time $X_{0}$:

$f_{\text{gate}}(x)\stackrel{{\scriptstyle\text{def}}}{{=}}(f_{\text{max}}*f_{0})(x)=\int\limits_{-\infty}^{\infty}f_{\text{max}}(x^{\prime})f_{0}(x-x^{\prime})dx^{\prime}.$
(12)

This convolution formula is general and valid for any PDFs, which do not have
to be Gaussian. On the other hand, this is the bottleneck of the gate-level
analysis: expression (12) can be calculated numerically (requiring
computational resources and time). Thus, if this formula is applied
numerically to every gate in a VLSI circuit, the analysis would not be
completed in acceptable time. For this reason, the challenge here is to find
an approach to approximate this formula and get some sort of an analytical
expression for convolution (12).

Figure 3: Illustration of delay propagation through a logic gate. At stage
(a), two signals arrive at the input of the gate. At stage (b), the
$\max$-operation is performed, which gives a skewed PDF. At the same time, the
gate has its own operation time described by some distribution (c). Thus, the
distribution of the gate delay (d) requires the convolution of the obtained
distribution (b) and given (c). This convolution results in a new
distribution, which clearly has a non-Gaussian form.

Let us look at the operation of a logic gate in detail. Figure 3 visualises
the step-by-step transformations of all the PDFs involved into the operation
of a logic gates. For the purpose of illustration, we consider that the
arrival and operation time distributions (i.e. the distributions of $X_{1}$,
$X_{2}$, and $X_{0}$) are Gaussian. These distributions do not have to be
Gaussian, but we want to show that even if _assumed_ Gaussian, the convolution
(12) will result into a non-Gaussian distribution, as clearly visualised in
that Figure.

Because the PDF for the output of a gate has obviously a non-Gaussian form and
does not imply any analytical solutions, the following ideas have circulated
within the Timing Analysis community. The traditional and overall accepted
approach is to approximate the $\max$ operation [3, 4]. We have thoroughly
discussed different approaches in the review Section III, where we have
grouped the block–based approaches into two sub-classes, with phenomenological
and parameterised delay models. Azuma et al. [15] have proposed to model the
PDF of the $\max$ with a GMM, which makes it the closest one to the present
study. However, the authors did not address the convolution in (12) and did
not explain how such a GMM decomposition should be obtained from an arbitrary
non-Gaussian distribution.

In this paper, firstly we derive the exact analytical formula for the gate
delay problem (12) that gives the PDF of (11). The assumption we make is that
the RVs $X_{i}$ ($i=0,1,2$) in (11) are Gaussian ones. Secondly, we propose to
handle non-Gaussian distributions as _linear combinations_ of functions of the
shape (1), as schematically shown in Figure 4. Such _basis_ functions are
called RBFs, and this technique is known as GMM and is widely used in signal
processing [60] and machine learning [61, 28]. We shall discuss Gaussian
Mixture Models in the next Chapter. Since the RBFs are Gaussian-like, our
approach allows one to traverse a delay’s PDF through a timing graph without
any loss of information by utilising the exact solution to (12). Note that a
common source for such a loss is due to approximation of non-Gaussian PDFs
with Gaussian distributions. Algorithms for decomposition a general non-
Gaussian distribution into a GMM and further its traversal will be presented
in Section VII. Now, we will focus on the analytical solution to the gate
delay problem (12) and theoretical foundations of our approach.

Figure 4: High-level illustration of the decomposition of a non-Gaussian
function and its representation by means of a sum of Radial Basis Functions.

##  V Distribution of the Gate Delay

###  V.1 Maximum of Two Correlated Gaussians

Before we move to the analysis, let us make a few useful remarks.

Firstly, we shall represent symmetric forms like the PDF for the maximum of
two correlated Gaussian RVs (6) as follows:

$g(x,\rho)=\sum\limits_{\begin{subarray}{c}i,j=1,2\\\ i\neq
j\end{subarray}}\phi_{i}(x)\Phi\left[\frac{1}{\sqrt{1-\rho^{2}}}\left(\frac{x-\mu_{j}}{\sigma_{j}}-\rho\frac{x-\mu_{i}}{\sigma_{i}}\right)\right].$
(13)

Figure 5: Comparison of different expressions for the PDF of the maximum of
two Gaussian RVs, $X_{1}\sim\mathcal{N}(1,0.5)$ and
$X_{2}\sim\mathcal{N}(3,3)$ for different correlation coefficients $\rho$. The
blue line shows the exact expression (13), the orange line shows the linear-
in-$\rho$ expression (14), and the dashed black line corresponds to the
independent case (5). For the purpose of illustration, $X_{1}$ and $X_{2}$ in
this example are two Gaussian RVs with different $\mu$ and $\sigma$.

Note that for $\rho=0$, we immediately obtain the PDF of the maximum of two
_independent_ Gaussian RVs, given by formula (5): $g(x,0)\equiv g(x)$.

Secondly, the key requirement for the algorithms discussed further is
_linearity_ , as we want to represent delay distributions via sums of RBFs.
For this purpose, we consider _weak correlations_ between RVs. In the case of
a small correlation coefficient, $|\rho|\rightarrow 0$, the function
$\Phi[\cdot]$ reads

$\Phi\left[\frac{1}{\sqrt{1-\rho^{2}}}\left(\frac{x-\mu_{j}}{\sigma_{j}}-\rho\frac{x-\mu_{i}}{\sigma_{i}}\right)\right]=\Phi\left(\frac{x-\mu_{j}}{\sigma_{j}}\right)-\phi_{j}(x)\frac{x-\mu_{i}}{\sigma_{i}}\sigma_{j}\rho+\ldots,$

and the PDF (13) can be approximated in terms of $\rho$ linearly:

$g(x,\rho)\approx g(x)+\rho\cdot\delta f_{\text{max}}(x),$ (14a) where $\delta
f_{\text{max}}(x)=-\phi_{1}(x)\phi_{2}(x)\left[\frac{x-\mu_{1}}{\sigma_{1}}\sigma_{2}+\frac{x-\mu_{2}}{\sigma_{2}}\sigma_{1}\right].$
(14b)

The formulas (13) and (14) are plotted in Figure 5 for
$X_{1}\sim\mathcal{N}(1,0.5)$ and $X_{2}\sim\mathcal{N}(3,3)$ and for
different correlations coefficients $\rho$. As a reference, the PDF of the
maximum or these RVs given by (5), when they are independent, i.e. neglecting
the correlations, is also shown. As one can see, the linear approximation (14)
performs well for a small correlation coefficient, however, it breaks down for
large $\rho$ as expected.

In the next Section, we derive the PDF for the gate delay distribution (12)
assuming initial delays to be Gaussian ones and arrival time delays to be
correlated. Then, we will relax this requirement and derive an approximation
for weak correlations. We also present the first four moments of the
distribution (12).

###  V.2 Derivation of the Exact Formula

Let us consider three Gaussian RVs, $X_{1}$, $X_{2}$ and $X_{0}$, and an RV
$\eta$ such that $\eta=\max(X_{1},X_{2})+X_{0}$. We shall assume that RVs
$X_{1}$ and $X_{2}$ are bounded with the correlation coefficient $\rho$, thus,
the PDF of $\max(X_{1},X_{2})$ is given by (13), and the PDF $\phi_{0}(x)$ of
the RV $X_{0}$ is given by (1). The variables $X_{1}$ and $X_{2}$ have the
meaning of the arrival times’ delays with mean values $\mu_{1}$, $\mu_{2}$ and
standard deviations $\sigma_{1}$, $\sigma_{2}$ correspondingly, i.e.
$X_{1}\sim\mathcal{N}(\mu_{1},\sigma_{1})$ and
$X_{2}\sim\mathcal{N}(\mu_{2},\sigma_{2})$, while
$X_{0}\sim\mathcal{N}(\mu_{0},\sigma_{0})$ . Having said that, we can write
(12) explicitly implying the Gaussian nature of its components $X_{i}$:

$f_{\text{gate}}(x)=\int\limits_{-\infty}^{\infty}g(x^{\prime},\rho)\phi_{0}(x-x^{\prime})dx^{\prime}.$
(15)

Substituting (13) into the latter equation, we obtain a sum of two integrals
of the form:

$I_{ij}^{(0)}=\int\limits_{-\infty}^{\infty}\phi_{i}(x^{\prime})\phi_{0}(x-x^{\prime})\times\Phi\left[\frac{1}{\sqrt{1-\rho^{2}}}\left(\frac{x^{\prime}-\mu_{j}}{\sigma_{j}}-\rho\frac{x^{\prime}-\mu_{i}}{\sigma_{i}}\right)\right]dx^{\prime}.$

Recalling the exact form of the PDF $\phi(x)$ and the error function
representation of the CDF $\Phi(x)$, formulas (1) and (2) respectively, and
introducing the intermediate variables,

$\displaystyle
a=\frac{1}{2}\frac{\sigma_{0}^{2}+\sigma_{i}^{2}}{\sigma_{0}^{2}\sigma_{i}^{2}},\quad
b=\frac{\sigma_{0}^{2}\mu_{i}+\sigma_{i}^{2}(x-\mu_{0})}{\sigma_{0}^{2}\sigma_{i}^{2}},\quad
c=-\frac{1}{2}\frac{\sigma_{0}^{2}\mu_{i}^{2}+\sigma_{i}^{2}(x-\mu_{0})^{2}}{\sigma_{0}^{2}\sigma_{i}^{2}},$

we arrive to the following expression:

$I_{ij}^{(0)}=\frac{\sqrt{2\pi}}{2}\sigma_{j}A_{0}A_{1}A_{2}\exp\left(\frac{b^{2}}{4a}+c\right)\\\
\times\int\limits_{-\infty}^{\infty}\exp\left[-a\left(x^{\prime}-\frac{b}{2a}\right)^{2}\right]\\\
\times\left\\{1+\operatorname{erf}\left[\frac{1}{\sqrt{2(1-\rho^{2})}}\left(\frac{x^{\prime}-\mu_{j}}{\sigma_{j}}-\rho\frac{x^{\prime}-\mu_{i}}{\sigma_{i}}\right)\right]\right\\}dx^{\prime}.$
(16)

Here, $A_{i}=1/\sqrt{2\pi}\sigma_{i}$ ($i=0,1,2$) are the normalisation
factors corresponding to the distributions $\phi_{i}(x)$. We present these
factors explicitly in order to investigate the most general case of Gaussian-
like distributions, which are not necessarily unit-normalised. The first
integral in (16) simply gives:

$\int\limits_{-\infty}^{\infty}e^{-a\left(x^{\prime}-\frac{b}{2a}\right)^{2}}dx^{\prime}=\sqrt{\frac{\pi}{a}}.$

Figure 6: PDF of the gate delay for a case of $\rho=0.5$. It is taken
$X_{1}\sim\mathcal{N}(1,0.75)$, $X_{2}\sim\mathcal{N}(2,3)$ and
$X_{0}\sim\mathcal{N}(0,1)$. Then, each of the parameters is varied as
specified keeping the rest ones constant.

The second integral in (16) can be brought to the form given by formula
(LABEL:eq:int_gaussian_with_erf) of Appendix LABEL:app:integrals. Indeed,
using a new variable of integration, $t=x^{\prime}-b/2a$, and introducing

$\displaystyle k$
$\displaystyle=\frac{1}{\sqrt{1-\rho^{2}}}\left(\frac{1}{\sigma_{j}}-\frac{\rho}{\sigma_{i}}\right),$
(17) $\displaystyle p$
$\displaystyle=\frac{\sigma_{0}^{2}\mu_{i}+\sigma_{i}^{2}(x-\mu_{0})}{\tilde{\sigma}_{i}^{2}}k-\frac{1}{\sqrt{1-\rho^{2}}}\left(\frac{\mu_{j}}{\sigma_{j}}-\rho\frac{\mu_{i}}{\sigma_{i}}\right),$

for the second integral, we have

$\int\limits_{-\infty}^{\infty}e^{-at^{2}}\operatorname{erf}\left(\frac{kt+p}{\sqrt{2}}\right)dt=\sqrt{\frac{\pi}{a}}\cdot\operatorname{erf}\left(\frac{p}{\sqrt{2}}\sqrt{\frac{a}{a+k^{2}/2}}\right).$

The integration procedure in the formula above is not trivial and discussed in
Appendix LABEL:app:integrals mentioned earlier.

Now we can write the formula for the distribution of the logic gate delay.
Noting that

$\sqrt{\frac{\pi}{a}}=\sqrt{2\pi}\frac{\sigma_{0}\sigma_{i}}{\sqrt{\sigma_{0}^{2}+\sigma_{i}^{2}}},\quad\frac{b^{2}}{4a}+c=-\frac{1}{2}\frac{(x-\mu_{0}-\mu_{i})^{2}}{\sigma_{0}^{2}+\sigma_{i}^{2}},$

we change to the initial variable of integration instead of the variables $k$
and $p$, and we return to the $\Phi$-representation. After simple algebraic
manipulations, we finally obtain:

$f_{\text{gate}}(x,\rho)=C\sum\limits_{\begin{subarray}{c}i,j=1,2\\\ i\neq
j\end{subarray}}\frac{1}{\tilde{\sigma}_{i}}\varphi\left(\frac{x-\mu_{0}-\mu_{i}}{\tilde{\sigma}_{i}}\right)\\\
\times\Phi\left[\frac{1}{\sqrt{1+\kappa_{ij}(\rho)^{2}}}\left(\frac{\sigma_{i}^{2}(x-\mu_{0})+\sigma_{0}^{2}\mu_{i}}{\tilde{\sigma}_{i}\sigma_{i}\sigma_{0}}\kappa_{ij}(\rho)-\frac{\sigma_{i}\mu_{j}-\rho\sigma_{j}\mu_{i}}{\sigma_{i}\sigma_{j}\sqrt{1-\rho^{2}}}\right)\right],$
(18a) where $C=2\pi A_{0}A_{1}A_{2}\sigma_{0}\sigma_{1}\sigma_{2}$ and
$\displaystyle\kappa_{ij}(\rho)$
$\displaystyle=\frac{\sigma_{i}-\rho\sigma_{j}}{\sqrt{1-\rho^{2}}}\cdot\frac{\sigma_{0}}{\sigma_{j}\tilde{\sigma}_{i}}=\frac{\kappa_{ij}}{\sqrt{1-\rho^{2}}}\left(1-\rho\frac{\sigma_{j}}{\sigma_{i}}\right),$
(18b) $\displaystyle\kappa_{ij}$
$\displaystyle=\frac{\sigma_{i}\sigma_{0}}{\sigma_{j}\tilde{\sigma}_{i}},\quad\tilde{\sigma}_{i}=\sqrt{\sigma_{0}^{2}+\sigma_{i}^{2}}\,.$
Separating $\kappa_{ij}$ explicitly in (18b), we can write the formula for the
PDF of the logic gate delay as follows:
$f_{\text{gate}}(x,\rho)=\,\,C\sum\limits_{\begin{subarray}{c}i,j=1,2\\\ i\neq
j\end{subarray}}\frac{1}{\tilde{\sigma}_{i}}\varphi\left(\frac{x-\mu_{0}-\mu_{i}}{\tilde{\sigma}_{i}}\right)\\\
\times\Phi\left\\{\left(\frac{1-\rho^{2}}{\left(1-\rho\sigma_{j}/\sigma_{i}\right)^{2}}+\kappa_{ij}^{2}\right)^{-\frac{1}{2}}\right.\\\
\times\left.\left[\eta_{ij}(x)\cdot\mathrm{sign}\left(1-\rho\frac{\sigma_{j}}{\sigma_{i}}\right)-\frac{\mu_{j}}{\sigma_{j}}\frac{1-\rho\sigma_{j}\mu_{i}/\sigma_{i}\mu_{j}}{|1-\rho\sigma_{j}/\sigma_{i}|}\right]\right\\}.$
(18c) For convenience in the formula above, we introduced the following
function:
$\eta_{ij}(x)=\frac{\sigma_{i}^{2}(x-\mu_{0})+\sigma_{0}^{2}\mu_{i}}{\tilde{\sigma}_{i}^{2}\sigma_{j}},$
(18d) and the $\mathrm{sign}$-function is defined in a conventional manner:
$\mathrm{sign}\left(1-\rho\frac{\sigma_{j}}{\sigma_{i}}\right)=\begin{cases}1,&1>\rho\sigma_{j}/\sigma_{i};\\\
0,&1=\rho\sigma_{j}/\sigma_{i};\\\
-1,&1<\rho\sigma_{j}/\sigma_{i}.\end{cases}$ (18e)

Let us discuss the obtained result. The formulas (18) give an _exact_ solution
to the gate delay problem (12) assuming that arrival times and gate operation
time (interconnect delay, etc.) are given by Gaussian distributions. The
resulting distribution (18) has a characteristic non-Gaussian form and is
fully determined by the parameters of the individual distributions of the
signals, i.e. by mean values $\mu_{i}$ and stds $\sigma_{i}$ of the arrival
times and gate operation time. The functions in (18), namely $\varphi(.)$ and
$\Phi(.)$, are the Gaussian kernel and Gaussian CDF given by equations (1) and
(2) respectively. These functions are well implemented in standard math
environments, which makes computation of (18) efficient.

The representations of the gate delay (18a) and (18c) are equivalent, and the
PDF is shown in Figure 6. However, from formula (18c) we can see that the
shape of a gate delay significantly simplifies if the input arrival delays (i)
have equal standard deviations, $\sigma_{1}=\sigma_{2}$, and (ii) are fully
correlated, $\rho=1$. If, in addition, the arrival delays have equal mean
values, $\mu_{1}=\mu_{2}$, the solution simplifies to a trivial case of the
convolution of two Gaussian distributions.

Figure 7: A family of curves for the PDF of the gate delay for different
values of the correlation coefficient $\rho$. As in Figure 6, the random
variables are taken as $X_{1}\sim\mathcal{N}(1,0.75)$,
$X_{2}\sim\mathcal{N}(2,3)$, and $X_{0}\sim\mathcal{N}(0,1)$.

###  V.3 Limit Cases of the Correlation Coefficient

Let us consider several cases with different strengths of the correlation
coefficient $\rho$. The family of curves for various values of $\rho$ is shown
in Figure 7. The ‘simplest’ result is achieved for independent arrival delays.
Hence, for $\rho=0$ the parameter

$\kappa_{ij}(\rho=0)=\kappa_{ij}=\frac{\sigma_{i}\sigma_{0}}{\sigma_{j}\tilde{\sigma}_{i}},$
(19)

and from (18a) one obtains:

$f_{\text{gate}}(x)=\frac{1}{\sqrt{2\pi}}\sum\limits_{\begin{subarray}{c}i,j=1,2\\\
i\neq
j\end{subarray}}\frac{1}{\tilde{\sigma}_{i}}\varphi\left(\frac{x-\mu_{0}-\mu_{i}}{\tilde{\sigma}_{i}}\right)\\\
\times\Phi\left[\frac{1}{\sqrt{1+\kappa_{ij}^{2}}}\left(\frac{\sigma_{i}^{2}(x-\mu_{0})+\sigma_{0}^{2}\mu_{i}}{\tilde{\sigma}_{i}^{2}\sigma_{j}}-\frac{\mu_{j}}{\sigma_{j}}\right)\right],$
(20)

where we also have written the value of the coefficient $C$ explicitly.

In the case when correlations are small, $\rho\rightarrow 0$, a small
parameter in the argument of the function $\Phi\\{\ldots\\}$ in (18a) appears,
and we can write:

$\Phi\left\\{\ldots\right\\}=\Phi\left[\frac{\eta_{ij}(x)-\mu_{j}/\sigma_{j}}{\sqrt{1+\kappa_{ij}^{2}}}\right]\\\
-\frac{1}{\sqrt{2\pi}(1+\kappa_{ij}^{2})^{3/2}}\frac{\sigma_{j}}{\sigma_{i}}\exp\left[\frac{(\eta_{ij}(x)-\mu_{j}/\sigma_{j})^{2}}{2(1+\kappa_{ij}^{2})}\right]\\\
\times\left(\eta_{ij}(x)-\frac{\mu_{j}}{\sigma_{j}}-\frac{(\mu_{i}-\mu_{j})(1+\kappa_{ij}^{2})}{\sigma_{j}}\right)\rho+\ldots$
(21)

Therefore, for the PDF (18a)–(18c) in the weak correlation case, we have:

$f_{\text{gate}}(x,\rho)\approx f_{\text{gate}}(x)+\rho\cdot\delta
f_{\text{gate}}(x),$ (22a) where $f_{\text{gate}}(x)$ is the PDF (20) for
independent arrival time delays, and the correction $\delta f_{\text{max}}$
reads: $\delta
f_{\text{gate}}(x)=-\,\,C\sum\limits_{\begin{subarray}{c}i,j=1,2\\\ i\neq
j\end{subarray}}\frac{1}{\tilde{\sigma}_{i}}\varphi\left(\frac{x-\mu_{0}-\mu_{i}}{\tilde{\sigma}_{i}}\right)\\\
\times\frac{1}{\sqrt{2\pi}(1+\kappa_{ij}^{2})^{3/2}}\frac{\sigma_{j}}{\sigma_{i}}\exp\left[\frac{(\eta_{ij}(x)-\mu_{j}/\sigma_{j})^{2}}{2(1+\kappa_{ij}^{2})}\right]\\\
\times\left(\eta_{ij}(x)-\frac{\mu_{j}}{\sigma_{j}}-\frac{(\mu_{i}-\mu_{j})(1+\kappa_{ij}^{2})}{\sigma_{j}}\right).$
(22b)

The obtained asymptotic result is _linear_ in terms of the correlation
coefficient $\rho$, which allows one to use it for the purposes of the
algorithm discussed in Section VII, where non-Gaussian delays are decomposed
into mixtures of RBFs. Here, let us examine the validity of this expansion for
different cases. We have already seen, that similar to (22) linear expression
(14) for the PDF of the maximum of two Gaussian RVs failed in strong
correlation regime (Figure 5). Therefore, let us discuss the issue of _strong_
correlations for the approximation $f_{\text{gate}}(x,\rho)$ below.

Figure 8: Comparison of the PDF of the gate delay obtained from different
expressions (the case of uncorrelated RVs, approximation and exact expression)
for four different correlation coefficients $\rho$. Gate operation time here
is $X_{0}\sim\mathcal{N}(3,0.5)$.

Figures 8–9 show a comparison of the exact result (18) with the linear in
$\rho$ approximation (22) for different values of the correlation coefficient
$\rho$. The random variables $X_{1,2}$ are taken with Gaussian distributions:
$X_{1}\sim\mathcal{N}(1,0.5)$ and $X_{2}\sim\mathcal{N}(3,1.9)$. In Figure 8,
the RV $X_{0}\sim\mathcal{N}(3,0.5)$, and we can see that the linear
approximation also fails when the correlations between the inputs are strong.
At the same time, Figure 9 shows the case with $X_{0}\sim\mathcal{N}(3,1)$,
i.e. with the twice larger standard deviation, and the linear approximation
(22) gives perfect agreement with the exact solution.

This should not be surprising, if one recalls that the gate delay PDF is the
_convolution_ of $g(x,\rho)$ with the Gaussian PDF $\phi_{0}(x)$, which
describes the delay due to the gate operation time, interconnect delay, _etc_.
Thus, when the variance of the latter is high, its contribution in the
resulting delay (15) becomes dominating, and inaccuracies in the approximation
of $g(x,\rho)$ vanish.

Figure 9: Comparison of the PDF of the gate delay obtained from different
expressions (the case of uncorrelated RVs, approximation and exact expression)
for four different correlation coefficients $\rho$. Gate operation time here
is $X_{0}\sim\mathcal{N}(3,1)$.

##  VI Moments of the Gate Delay Distribution

In this Section, we consider the first moments of the distribution (18a). In
general, the moments are given by the following expressions:

$\mathbb{E}[X^{n}]=\int\limits_{-\infty}^{\infty}x^{n}f_{\text{gate}}(x,\rho)dx.$
(23)

Below, we discuss the first four moments: the mean, standard deviation,
skewness and kurtosis of the gate delay PDF. The meaning of these moments is
explained in Section LABEL:chap2:sec:basic_prob.

###  VI.1 Mean

The mean of the gate delay is given by

$\mu_{\text{gate}}\equiv\mathbb{E}[X]=\int\limits_{-\infty}^{\infty}xf_{\text{gate}}(x,\rho)dx.$
(24)

The substitution of (18a) into (24) leads to the calculations, similar to
those performed above. Therefore, in this Section we only present the final
results for the moments. The mean value of the gate delay reads:

$\mu_{\text{gate}}=C\sum\limits_{\begin{subarray}{c}i,j=1,2\\\ i\neq
j\end{subarray}}\left\\{\sqrt{2\pi}(\mu_{0}+\mu_{i})\Phi\left[\frac{1}{\sqrt{1-\rho^{2}}}\frac{\mu_{i}-\mu_{j}}{\sigma_{j}\varkappa_{ij}}\right]\right.\\\
+\left.\frac{\sigma_{i}\tilde{\sigma}_{i}}{\sigma_{0}}\frac{\kappa_{ij}(\rho)}{\varkappa_{ij}}\exp\left[-\frac{1}{2}\frac{1}{1-\rho^{2}}\frac{(\mu_{i}-\mu_{j})^{2}}{\sigma_{j}^{2}\varkappa_{ij}^{2}}\right]\right\\},$
(25)

where

$\varkappa_{ij}=\sqrt{1+\kappa_{ij}(\rho)^{2}\left(1+\frac{\sigma_{i}^{2}}{\sigma_{0}^{2}}\right)}.$
(26)

Figure 10 shows the dependencies of the mean value $\mu_{\text{gate}}$ on the
mean values of input delays and the standard deviations of input delays. The
gate operation time is taken as a standard Gaussian
$X_{0}\sim\mathcal{N}(0,1)$, and the correlation coefficient is $\rho=0.5$.

###  VI.2 Standard Deviation

The evaluation of the standard deviation requires the knowledge of the second
central moment:

$\sigma_{\text{gate}}\equiv\mathbf{D}X\equiv\sqrt{\langle X^{2}\rangle-\langle
X\rangle^{2}}=\sqrt{\mathbb{E}(X^{2})-\mathbb{E}(X)^{2}}.$ (27)

Performing the integration in (23) for $n=2$, after some algebra, for the
second central moment one obtains:

$\mathbb{E}[X^{2}]=C\sum\limits_{\begin{subarray}{c}i,j=1,2\\\ i\neq
j\end{subarray}}\left\\{\sqrt{2\pi}\left[(\mu_{0}+\mu_{i})^{2}+\tilde{\sigma}_{i}^{2}\right]\Phi\left[\frac{y(\ldots)}{\varkappa_{ij}}\right]\right.\\\
+2\frac{\sigma_{i}\tilde{\sigma}_{i}^{2}}{\sigma_{0}}\frac{\kappa_{ij}(\rho)}{\varkappa_{ij}}\left(\frac{\mu_{0}+\mu_{i}}{\tilde{\sigma}_{i}}-\frac{1}{2}\frac{\sigma_{i}}{\sigma_{0}}\frac{\kappa_{ij}(\rho)}{\varkappa_{ij}^{2}}y(\ldots)\right)\\\
\times\left.\exp\left[-\frac{1}{2}\frac{y(\ldots)^{2}}{\varkappa_{ij}^{2}}\right]\right\\},$
(28)

where

$y(\ldots)=\frac{\tilde{\sigma}_{i}\mu_{i}}{\sigma_{i}\sigma_{0}}\kappa_{ij}(\rho)-\frac{\mu_{j}}{\sigma_{j}}\frac{1-\rho\sigma_{j}\mu_{i}/\sigma_{i}\mu_{j}}{\sqrt{1-\rho^{2}}}\,.$
(29)

Now, having both first and second central moments, $\mathbb{E}[X]$ and
$\mathbb{E}[X^{2}]$, we can calculate the standard deviation
$\sigma_{\text{gate}}$ of the gate delay. This quantity is presented in Figure
11 for the same values of the parameters as in the Figure with the mean value
$\mu_{\text{gate}}$.

Figure 10: Mean value of the gate delay as a function of the input delay mean
values and standard deviations for $\rho=0.5$ and $X_{0}\sim\mathcal{N}(0,1)$.
Figure 11: Gate delay standard deviation as a function of $\mu_{1,2}$ and
$\sigma_{1,2}$. The parameters are the same as in Figure 10.

###  VI.3 Higher Moments

The higher moments can be obtained in the same manner. We do not present the
exact expressions due to their cumbersomeness. However, gate delay skewness

$\gamma_{\text{gate}}\equiv\operatorname{E}\left[\left(\frac{X-\mu}{\sigma}\right)^{3}\right]=\frac{\operatorname{E}[X^{3}]-3\mu_{\text{gate}}\sigma_{\text{gate}}^{2}-\mu_{\text{gate}}^{3}}{\sigma_{\text{gate}}^{3}},$
(30)

and the kurtosis

$\kappa_{\text{gate}}\equiv\operatorname{E}\left[\left(\frac{X-\mu}{\sigma}\right)^{4}\right]=\frac{\operatorname{E}[(X-\mu)^{4}]}{\sigma^{4}},$
(31)

are shown in Figures 12–13 correspondingly.

Figure 12: Gate delay skewness as the function of input delays’ mean values
and standard deviations for $\rho=0.5$ and $X_{0}\sim\mathcal{N}(0,1)$. Figure
13: Gate delay kurtosis as the function of input delays’ mean values and
standard deviations for $\rho=0.5$ and $X_{0}\sim\mathcal{N}(0,1)$.

##  VII Delay Propagation Algorithm

###  VII.1 Brief Theoretical Summary

Consider a logic gate with two inputs, $A$ and $B$, and suppose that the gate
operation time is distributed according to the normal law with a mean
$\mu_{0}$ and a variance $\sigma_{0}^{2}$, i.e. it is a Gaussian RV. Assume
now that the arrival times of both signals are also Gaussian RVs with means
and variances ($\mu_{1}$, $\sigma_{1}^{2}$) and ($\mu_{2}$, $\sigma_{2}^{2}$)
respectively. Even if the individual distribution in the example of Figure 3
are Gaussian, the application of formula (11) leads to a non-Gaussian
distribution of the delay at the gate output as shown in that figure. The
exact expression for the PDF of such an RV is

$f_{\text{gate}}(x)=\frac{1}{\sqrt{2\pi}}\sum\limits_{\begin{subarray}{c}i,j=1,2\\\
i\neq
j\end{subarray}}\frac{1}{\tilde{\sigma}_{i}}\varphi\left(\frac{x-\mu_{0}-\mu_{i}}{\tilde{\sigma}_{i}}\right)\Phi\left[\frac{1}{\sqrt{1+\kappa_{ij}^{2}}}y(\ldots)\right],$
(32)

where

$\kappa_{ij}=\frac{\sigma_{0}\sigma_{i}}{\sigma_{j}\tilde{\sigma}_{i}},\quad\tilde{\sigma}_{i}=\sqrt{\sigma_{0}^{2}+\sigma_{i}^{2}},\quad
y(x)=\frac{\sigma_{i}^{2}(x-\mu_{0})+\sigma_{0}^{2}\mu_{i}}{\tilde{\sigma}_{i}^{2}\sigma_{j}}-\frac{\mu_{j}}{\sigma_{j}}.$
(33)

While expression (32) does not take into account possible correlations between
the arrival signals, we are interested in _demonstrating_ how this exact
solution can speed up SSTA for a given graph keeping precision high. The full
expression that takes the correlations into account was derived above and is
given by equation (18). Formula (32) assumes all initial delays (arrival and
gate itself) to have Gaussian distributions. In principle, both the arrival
signal and gate delay do not have to be Gaussian. If they can be decomposed
into a linear superposition of Gaussian kernel functions, the PDF of the gate
output delay can be presented as a linear combination of expressions (32) due
to the linearity of the integration operation. This idea constitutes the core
of a delay propagation algorithm which we discuss in the next Sections.

###  VII.2 Algorithm

The algorithm GaKeDA for the calculation of the delay propagation through a
timing graph is shown in Figure 14 and outlined below. The key feature of the
algorithm is representing non-Gaussian distributions with a mixture of RBFs
that have Gaussian form, i.e. as a linear combination (weighted sum) of
Gaussian kernels $\varphi(x)$ given by (1). This allows one to use the result
(32), obtained above for Gaussian distributions. Hence the name for the
algorithm: Gaussian Kernel Density Estimation based Algorithm (GaKeDA). The
Algorithm returns a list with the parameters for the Gaussian mixtures for
each node of the graph $G$.

Data: graph $G$ and delay distributions for nodes and input signals

Result: a list with PDFs for all nodes in the RBF representation

1 /* We assume that the sequence, in which the nodes should be visited, is
known */

for  _node in $G$_ do

$f_{node}\leftarrow$ initialise;

if  _delay’s PDF not Gaussian_ then

$f_{component}\leftarrow$ utilise the exact formula (32);

$f_{node}\leftarrow f_{node}+f_{component}$;

end if

else

$f_{node}\leftarrow$ utilise the exact formula (32); /* delay is Gaussian */

end if

$f_{GMM}\leftarrow$ decompose $f_{node}$ into GMM ;

add GMM parameters to a list

end for

Figure 14: GaKeDA Algorithm.

Such models that use Gaussian kernels (1) as the shape functions are called
GMMs and are well known and studied in the literatureThere are numerous books
and monographs on mixture models, not only GMMs, that cover theoretical
foundations and applications, e.g. by Everitt & Hand [63], Titterington et al.
[64], McLachlan & Basford [65], Lindsay [66], Böhning [67], McLachlan & Peel
[68], Frühwirth–Schnatter [69], Mengersen et al. [70], and McNicholas [71]. A
recent review by McLachlan et al. [72] was published in 2019. Applications of
mixture models in Machine Learning are covered in a book by Murphy [28]. GMMs
are realised in various software packages such as scikit-learn [73].

In case the arrival signals’ delays and operation time of a gate have Gaussian
distribution, the formula (32) is applied straightforwardly (see the line $9$
of GaKeDA in Figure 14). When one of the delays has non-Gaussian distribution,
it should be first decomposed into a Gaussian mixture and, therefore, the
formula (32) is applied for every pair of the mixture components (see the
lines $5$ – $7$ of GaKeDA). From here, one can see the complexity of the
Algorithm: in worst case, when all the delays in $G$ are represented via
Gaussian mixtures, the complexity with respect to the total number $N$ of the
atomic operations is $\mathcal{O}(m^{3}N)$, where $m$ is a number of
components in the Gaussian mixture. However, the gate operation time often has
Gaussian distribution [74]. In this case, only the arrival times should be
represented by means of the mixtures and, thus, the complexity reduces to
$\mathcal{O}(m^{2}N)$. Such complexity is comparable with that of the known
approaches within the block–based SSTA, which we have discussed in Section III
(see Tables 1 and 2)

The high-level diagram of the GaKeDA flow is shown in Figure 15. The algorithm
relies on the decomposition procedure (see the line $10$ of GaKeDA) to find
the parameters of GMM representation. The GMM parameters identification
problem is a non-convex optimisation problem, thus, GaKeDA may be inefficient.
In the next section, we propose a Gaussian comb model for GMM, which allows
efficient GMM parameters identification.

Figure 15: High-level diagram of the GaKeDA flow.

###  VII.3 Optimisation

The exact function (32) can be written as $f_{\text{rbf}}(x)$, a sum of RBFs,
and each of these RBFs has a Gaussian-like shape. In other words, it can be
decomposed into a Gaussian mixture [75, 64, 68, 76]. The decomposition
procedure is equivalent to fitting the actual PDF with a sum of RBFs, which
brings us to an optimisation problem. In this study, we discuss the
minimisation of the sum of squares of the residuals

$\min\sum_{i}|f_{\text{rbf}}(x_{i})-y_{i}|^{2},$ (34)

subject to constraints (specified below). Here $y_{i}$ are the data points
corresponding to the actual function $f_{\text{gate}}(x)$ that we want to fit.

Depending on the form of the RBFs, the minimisation of (34) can vary
significantly, e.g. an approximate function $f_{\text{rbf}}(x)$ can be either
linearly or non-linearly dependent on the fitting parameters. Different
approaches to this problem were discussed by Mishagli & Blokhina [7].

The optimisation problem should be formulated as a LP problem in order to
avoid a numerical approximate optimisation procedure and to speed up the
algorithm [77]. Instead of the $2-$norm $||\ldots||_{2}$, consider the
infinite norm

$||\mathbb{M}\mathbf{x}-\mathbf{y}||_{\infty}=\max(|r_{1}|,|r_{2}|,\ldots,|r_{n}|),$
(35)

where $|r_{i}|=|M_{i}^{T}\mathbf{x}-y_{i}|$ ($i=1,\ldots,n$) are residuals.
Denoting $\max(\ldots)$ by $t$, we have the following LP problem:

$\begin{array}[]{ll}\min&t\\\
\text{s.t.}&-t\mathbf{1}\preceq\mathbb{M}\mathbf{x}-\mathbf{y}\preceq
t\mathbf{1},\quad\mathbf{x}\succeq\mathbf{0}.\end{array}$ (36)

Thus, for an LP problem of the form

$\begin{array}[]{ll}\min&\mathbf{c}^{T}\mathbf{x}\\\
\text{s.t.}&\mathbb{A}\mathbf{x}\preceq\mathbf{b},\quad\mathbf{x}\succeq\mathbf{0},\end{array}$
(37)

the vectors $\mathbf{c}$, $\mathbf{x}$ and $\mathbf{b}$ and the matrix
$\mathbf{A}$ are as follows. Having $n$ data points ($x_{i},y_{i}$) and
considering $m$ RBFs, the vectors $\mathbf{c}$ and $\mathbf{x}$ are

$\mathbf{c}^{T}=(\underbrace{0,\ldots,0}_{m\text{
times}},1)\quad\mathbf{x}^{T}=(w_{1},w_{2},\ldots,w_{m},t),$ (38)

where $w_{i}$ ($i=1,\dots,m$) are unknown weights.

The matrix $\mathbb{A}$ of size $(2n,m+1)$ reads

$\mathbb{A}=\begin{pmatrix}\mathbb{M}&-\mathbf{1}\\\
-\mathbb{M}&-\mathbf{1}\end{pmatrix},$ (39)

where

$\mathbb{M}=\begin{pmatrix}\varphi_{1}(x_{1})&\varphi_{2}(x_{1})&\ldots&\varphi_{m}(x_{1})\\\
\varphi_{1}(x_{2})&\varphi_{2}(x_{2})&\ldots&\varphi_{m}(x_{2})\\\
\vdots&\vdots&\ddots&\vdots\\\
\varphi_{1}(x_{n})&\varphi_{2}(x_{n})&\ldots&\varphi_{m}(x_{n})\end{pmatrix}.$
(40)

The vector $\mathbf{b}^{T}=(\mathbf{y},-\mathbf{y})$ is of length $2n$.

##  VIII Discussion and Conclusions

In this paper, original research contributions in the context of block–based
SSTA were presented. The SSTA problem of a logic circuit delay computation was
considered at the gate level. A fundamental (mathematical) operation that
describes the operation of a logic gate, i.e. the _convolution_ (12), was
studied. To summarise, the main results of the paper are as follows.

  * •

The exact analytical solution to (12) has been given in form of (18). The only
assumption made was that the arrival times delays and gate operation times
were normally distributed. The expression for the gate delay PDF (18) is a
function of input distributions parameters, i.e. mean values, stds, and the
correlation coefficient of the arrival times; it has no unknown fitting
parameters. It should be noted that even with the assumption of normally
distributed input parameters, such a solution was not known in the field.

  * •

The obtained expression for the gate delay has the advantage to the known
results, such as by Azuma et al. [15], as the convolution is also taken into
account, thus, the gate delay problem is fully analysed.

  * •

The obtained solution has been studied in both zero and weak correlations
cases. In the latter case, the PDF of the gate delay is linear in the
correlation coefficient $\rho$. It has also been shown that, depending on the
standard deviation of gates operation times, this approximation can be
applicable in the strong correlation regime as well.

  * •

For the derived distribution, the first and second moments, mean and standard
deviation, as well as the higher moments, skewness and kurtosis, have been
obtained.

  * •

The theoretical foundations for algorithms that rely on _decomposition_ of
non-Gaussian PDF into a superposition of RBFs have been given.

  * •

The exact formula for an output logic gate delay, the convolution of
$\max(X_{1},$ $X_{2})$ and $X_{0}$ for Gaussian RVs $X_{i}$ ($i=0,1,2$),
allows one to build a closed-loop algorithm of forward traversal of a delay
through a timing graph $G$. The requirement for this is that non-Gaussian PDFs
of delays are presented via Gaussian mixtures, sums of RBFs of Gaussian form.
The decomposition is equivalent to solving the minimisation problem (34).

  * •

The low complexity of the Gaussian comb allows construction of efficient SSTA
algorithms. Thus, Azumaet al. [15] also proposed to use GMMs but did not
propose a procedure for GMM preparation; Cheng et al. [33] used Fourier Series
decomposition for evaluation of PDFs, which required exhaustive numerical
computations. This is not the case for the proposed approach in this study.

## References

  * Sapatnekar [2004] S. Sapatnekar,  _Timing_ (Springer-Verlag, 2004).
  * Gerez [1998] S. H. Gerez,  _Algorithms for VLSI Design Automation_ (Wiley, 1998).
  * Blaauw  _et al._ [2008] 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. 4 (2008), 10.1109/TCAD.2007.907047.
  * Forzan and Pandini [2009] C. Forzan and D. Pandini, “Statistical static timing analysis: A survey,” Integration, the VLSI Journal 42, 409–435 (2009).
  * Lewis and Schmit [2019] D. Lewis and H. Schmit, “Spatial Timing Analysis With Exact Propagation of Delay and Application to FPGA Performance,” IEEE Trans. Comput.–Aided Des. Integr. Circuits Syst. 64, 2153–2166 (2019).
  * Freeley  _et al._ [2018] J. Freeley, D. Mishagli, T. Brazil, and E. Blokhina, “Statistical Simulations of Delay Propagation in Large Scale Circuits Using Graph Traversal and Kernel Function Decomposition,” in  _Proc. SMACD_ (2018).
  * Mishagli [2020] D. Mishagli, “RBF approximation for non-Gaussian delay propagation in VLSI: Source Code,” doi:10.5281/zenodo.374975 (2020).
  * Pillage and Rohrer [1990] L. Pillage and R. Rohrer, “Asymptotic waveform evaluation for timing analysis,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 9, 352–366 (1990).
  * Chiprout and Nakhla [1994] E. Chiprout and M. S. Nakhla,  _Asymptotic Waveform Evaluation_ (Springer US, Boston, MA, 1994).
  * Berkelaar [1997] M. R. C. M. Berkelaar, “Statistical delay calculation: A linear time method,” in  _Proceedings of TAU 97: ACM/IEEE International Workshop on Timing Issues in the Specification and Synthesis of Digital Systems_ (1997) pp. 15–24.
  * Clark [1961] C. E. Clark, “The Greatest of a Finite Set of Random Variables,” Oper. Res. 9, 145–162 (1961).
  * Jacobs and Berkelaar [2000] E. T. A. F. Jacobs and M. R. C. M. Berkelaar, “Gate sizing using a statistical delay model,” in  _Proceedings Design, Automation and Test in Europe Conference and Exhibition 2000_ (2000) pp. 283–290.
  * Tsukiyama, Tanaka, and Fukui [2000] S. Tsukiyama, M. Tanaka, and M. Fukui, “A new statistical static timing analyzer considering correlation between delays,” in  _Proceedings of TAU ’00: International Workshop on Timing Issues in the Specification and Synthesis of Digital Systems_ (2000) pp. 27–33.
  * Tsukiyama, Tanaka, and Fukui [2001] S. Tsukiyama, M. Tanaka, and M. Fukui, “A statistical static timing analysis considering correlations between delays,” in  _Proceedings of the ASP-DAC 2001_ (2001) pp. 353–358.
  * Azuma, Tsukiyama, and Fukui [2017] D. Azuma, S. Tsukiyama, and M. Fukui, “Approximating the maximum of Gaussians by a Gaussian mixture model for statistical designs,” in  _2017 European Conference on Circuit Theory and Design (ECCTD)_ (2017).
  * Chen  _et al._ [2015] J. Chen, S. Cotofana, S. Grandhi, C. Spagnol, and E. Popovici, “Inverse Gaussian distribution based timing analysis of Sub-threshold CMOS circuits,” Microelectronics Reliability 55, 2754–2761 (2015).
  * Jin  _et al._ [2022] L. Jin, W. Fu, H. Yan, and L. Shi, “A Statistical Cell Delay Model for Estimating the 3$\sigma$ Delay by Matching Kurtosis,” IEEE Transactions on Circuits and Systems II: Express Briefs 69, 2932–2936 (2022).
  * Chang and Sapatnekar [2003] H. Chang and S. S. Sapatnekar, “Statistical timing analysis considering spatial correlations using a single PERT-like traversal,” in  _Proc. ICCAD_ (2003) pp. 621–625.
  * Chang and Sapatnekar [2005] H. Chang and S. Sapatnekar, “Statistical timing analysis under spatial correlations,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 24, 1467–1482 (2005).
  * Chang  _et al._ [2005] H. Chang, V. Zolotov, S. Narayan, and C. Visweswariah, “Parameterized block-based statistical timing analysis with non-gaussian parameters, nonlinear delay functions,” in  _Proceedings of the 42nd Annual Design Automation Conference_, DAC ’05 (Association for Computing Machinery, New York, NY, USA, 2005) pp. 71–76.
  * Visweswariah  _et al._ [2004] C. Visweswariah, K. Ravindran, K. Kalafala, S. G. Walker, and S. Narayan, “First-Order Incremental Block-Based Statistical Timing Analysis,” in  _Proc. DAC_ (ACM, New York, NY, USA, 2004) pp. 331–336.
  * Visweswariah  _et al._ [2006] C. Visweswariah, K. Ravindran, K. Kalafala, S. G. Walker, S. Narayan, D. K. Beece, J. Piaget, N. Venkateswaran, and J. G. Hemmett, “First-Order Incremental Block-Based Statistical Timing Analysis,” IEEE Trans. Comput.–Aided Des. Integr. Circuits Syst. 25, 2170–2180 (2006).
  * Agarwal  _et al._ [2003] A. Agarwal, D. Blaauw, V. Zolotov, S. Sundareswaran, M. Zhao, K. Gala, and R. Panda, “Statistical delay computation considering spatial correlations,” in  _Proceedings of the ASP-DAC Asia and South Pacific Design Automation Conference, 2003._ (2003) pp. 271–276.
  * Singh and Sapatnekar [2006] J. Singh and S. Sapatnekar, “Statistical timing analysis with correlated non-Gaussian parameters using independent component analysis,” in  _2006 43rd ACM/IEEE Design Automation Conference_ (2006) pp. 155–160.
  * Singh and Sapatnekar [2008] J. Singh and S. S. Sapatnekar, “A Scalable Statistical Static Timing Analyzer Incorporating Correlated Non-Gaussian and Gaussian Parameter Variations,” IEEE Trans. Comput.–Aided Des. Integr. Circuits Syst. 27, 160–173 (2008).
  * Manduchi and Portilla [1999] R. Manduchi and J. Portilla, “Independent component analysis of textures,” in  _Proceedings of the Seventh IEEE International Conference on Computer Vision_, Vol. 2 (1999) pp. 1054–1060 vol.2.
  * Hyvärinen and Oja [2000] A. Hyvärinen and E. Oja, “Independent component analysis: Algorithms and applications,” Neural Networks 13, 411–430 (2000).
  * Murphy [2022] K. P. Murphy,  _Probabilistic Machine Learning: An Introduction_ , Adaptive Computation and Machine Learning Series (The MIT Press, Cambridge, Massachusetts, 2022).
  * Li  _et al._ [2004] X. Li, J. Le, P. Gopalakrishnan, and L. T. Pileggi, “Asymptotic probability extraction for non-normal distributions of circuit performance,” in  _Proc. ICCAD_ (San Jose, CA, USA, 2004) pp. 2–9.
  * Zhan  _et al._ [2005] Y. Zhan, A. J. Strojwas, X. Li, L. T. Pileggi, D. Newmark, and M. Sharma, “Correlation-aware statistical timing analysis with non-Gaussian delay distributions,” in  _Proc. DAC_ (2005) pp. 77–82.
  * Zhang  _et al._ [2005] L. Zhang, W. Chen, Y. Hu, J. A. Gubner, and C. C. P. Chen, “Correlation-preserved non-Gaussian statistical timing analysis with quadratic timing model,” in  _Proc. DAC_ (2005) pp. 83–88.
  * Khandelwal and Srivastava [2005] V. Khandelwal and A. Srivastava, “A general framework for accurate statistical timing analysis considering correlations,” in  _Proc. DAC_ (2005) pp. 89–94.
  * Cheng, Xiong, and He [2007] L. Cheng, J. Xiong, and L. He, “Non-Linear Statistical Static Timing Analysis for Non-Gaussian Variation Sources,” in  _2007 44th ACM/IEEE Design Automation Conference_ (2007) pp. 250–255.
  * Cheng  _et al._ [2012] L. Cheng, F. Gong, W. Xu, J. Xiong, L. He, and M. Sarrafzadeh, “Fourier Series Approximation for Max Operation in Non-Gaussian and Quadratic Statistical Static Timing Analysis,” IEEE Transactions on Very Large Scale Integration (VLSI) Systems 20, 1383–1391 (2012).
  * Cheng, Xiong, and He [2008] L. Cheng, J. Xiong, and L. He, “Non-Gaussian statistical timing analysis using second-order polynomial fitting,” in  _2008 Asia and South Pacific Design Automation Conference_ (2008) pp. 298–303.
  * Cheng, Xiong, and He [2009] L. Cheng, J. Xiong, and L. He, “Non-Gaussian Statistical Timing Analysis Using Second-Order Polynomial Fitting,” IEEE Trans. Comput.–Aided Des. Integr. Circuits Syst. 28, 130–140 (2009).
  * Vijaykumar and Vasudevan [2014] M. Vijaykumar and V. Vasudevan, “Statistical static timing analysis using a skew-normal canonical delay model,” in  _2014 Design, Automation Test in Europe Conference Exhibition (DATE)_ (2014).
  * Ramprasath, Vijaykumar, and Vasudevan [2016] S. Ramprasath, M. Vijaykumar, and V. Vasudevan, “A Skew-Normal Canonical Model for Statistical Static Timing Analysis,” IEEE Transactions on Very Large Scale Integration (VLSI) Systems 24, 2359–2368 (2016).
  * Shebaita, Le, and Ding [2018] A. Shebaita, J. Le, and L. Ding, “Generalized Moment Based Framework for Accurate and Efficient SSTA,” in  _Proc. TAU_ (2018) pp. 37–42.
  * Zhang, Shao, and Chen [2006] L. Zhang, J. Shao, and C. C.-P. Chen, “Non-Gaussian statistical parameter modeling for SSTA with confidence interval analysis,” in  _Proceedings of the 2006 International Symposium on Physical Design_, ISPD ’06 (Association for Computing Machinery, New York, NY, USA, 2006) pp. 33–38.
  * Chopra  _et al._ [2006] K. Chopra, B. Zhai, D. Blaauw, and D. Sylvester, “A New Statistical Max Operation for Propagating Skewness in Statistical Timing Analysis,” in  _2006 IEEE/ACM International Conference on Computer Aided Design_ (2006) pp. 237–243.
  * Azzalini [1985] A. Azzalini, “A Class of Distributions Which Includes the Normal Ones,” Scandinavian Journal of Statistics 12, 171–178 (1985).
  * Azzalini [2013] A. Azzalini,  _The Skew-Normal and Related Families_, Institute of Mathematical Statistics Monographs (Cambridge University Press, Cambridge, 2013).
  * Owen [1956] D. B. Owen, “Tables for Computing Bivariate Normal Probabilities,” The Annals of Mathematical Statistics 27, 1075–1090 (1956).
  * Agarwal, Chopra, and Blaauw [2005] A. Agarwal, K. Chopra, and D. Blaauw, “Statistical timing based optimization using gate sizing,” in  _Design, Automation and Test in Europe_ (2005) pp. 400–405 Vol. 1.
  * Jin  _et al._ [2019] L. Jin, W. Fu, Y. Zheng, and H. Yan, “A Precise Block-Based Statistical Timing Analysis with MAX Approximation Using Multivariate Adaptive Regression Splines,” in  _2019 IEEE 13th International Conference on ASIC (ASICON)_ (2019).
  * Zolotov  _et al._ [2007] V. Zolotov, J. Xiong, S. Abbaspour, D. J. Hathaway, and C. Visweswariah, “Compact modeling of variational waveforms,” in  _Proc ICCAD_ (2007) pp. 705–712.
  * Fatemi, Nazarian, and Pedram [2006] H. Fatemi, S. Nazarian, and M. Pedram, “Statistical logic cell delay analysis using a current-based model,” in  _Proc. DAC_ (SanFrancisco, California, USA, 2006) pp. 253–256.
  * Keller, Tseng, and Verghese [2004] I. Keller, K. Tseng, and N. Verghese, “A robust cell-level crosstalk delay change analysis,” in  _IEEE/ACM International Conference on Computer Aided Design, 2004. ICCAD-2004._ (2004) pp. 147–154.
  * Cao and Clark [2005] Y. Cao and L. Clark, “Mapping statistical process variations toward circuit performance variability: An analytical modeling approach,” in  _Proceedings. 42nd Design Automation Conference, 2005._ (2005) pp. 658–663.
  * Sakurai and Newton [1990] T. Sakurai and A. Newton, “Alpha-power law MOSFET model and its applications to CMOS inverter delay and other formulas,” IEEE Journal of Solid-State Circuits 25, 584–594 (1990).
  * Orshansky, Chen, and Hu [1999] M. Orshansky, J. Chen, and C. Hu, “Direct sampling methodology for statistical analysis of scaled CMOS technologies,” IEEE Transactions on Semiconductor Manufacturing 12, 403–408 (1999).
  * Shinkai  _et al._ [2006] K. Shinkai, M. Hashimoto, A. Kurokawa, and T. Onoye, “A Gate Delay Model Focusing on Current Fluctuation over Wide-Range of Process and Environmental Variability,” in  _2006 IEEE/ACM International Conference on Computer Aided Design_ (2006) pp. 47–53.
  * Shinkai, Hashimoto, and Onoye [2013] K.-i. Shinkai, M. Hashimoto, and T. Onoye, “A gate-delay model focusing on current fluctuation over wide range of process–voltage–temperature variations,” Integration 46, 345–358 (2013).
  * Sinha  _et al._ [2016] D. Sinha, V. Zolotov, S. K. Raghunathan, M. H. Wood, and K. Kalafala, “Practical statistical static timing analysis with current source models,” in  _Proc. DAC_ (2016) pp. 1–6.
  * Zolotov  _et al._ [2012] V. Zolotov, D. Sinha, J. Hemmett, E. Foreman, C. Visweswariah, J. Xiong, J. Leitzen, and N. Venkateswaran, “Timing analysis with nonseparable statistical and deterministic variations,” in  _DAC Design Automation Conference 2012_ (2012) pp. 1061–1066.
  * Cao  _et al._ [2019] P. Cao, Z. Liu, J. Guo, and J. Wu, “An Analytical Gate Delay Model in Near/Subthreshold Domain Considering Process Variation,” IEEE Access 7, 171515–171524 (2019).
  * Harris  _et al._ [2010] D. M. Harris, B. Keller, J. Karl, and S. Keller, “A transregional model for near-threshold circuits with application to minimum-energy operation,” in  _2010 International Conference on Microelectronics_ (2010) pp. 64–67.
  * Keller, Harris, and Martin [2014] S. Keller, D. M. Harris, and A. J. Martin, “A Compact Transregional Model for Digital CMOS Circuits Operating Near Threshold,” IEEE Transactions on Very Large Scale Integration (VLSI) Systems 22, 2041–2053 (2014).
  * Kay [1993] S. M. Kay,  _Fundamentals of Statistical Signal Processing: Estimation Theory_ (Prentice Hall, 1993).
  * Kelleher, Mac Namee, and D’Arcy [2015] J. D. Kelleher, B. Mac Namee, and A. D’Arcy,  _Fundamentals of Machine Learning for Predective Data Analytics_ (The MIT Press, 2015).
  * Note [1] The first known example of application of GMM for data modelling was done by Pearson in his 1894 paper[78].
  * Everitt and Hand [1981] B. Everitt and D. Hand,  _Finite Mixture Distributions_ (Chapman & Hall, New York, 1981).
  * Titterington, Smith, and Makov [1985] D. Titterington, A. Smith, and U. Makov,  _Statistical Analysis of Finite Mixture Distributions_ , Wiley Series in Probability and Mathematical Statistics: Applied Probability and Statistics (Wiley, 1985).
  * McLachlan and Basford [1987] G. J. McLachlan and K. E. Basford,  _Mixture Models: Inference and Applications to Clustering_ , Statistics: A Series of Textbooks and Monographs, Vol. 84 (Taylor & Francis, 1987).
  * Lindsay [1995] B. Lindsay,  _Mixture Models: Theory, Geometry, and Applications_ , Conference Board of the Mathematical Sciences: NSF-CBMS Regional Conference Series in Probability and Statistics (Institute of Mathematical Statistics, 1995).
  * Böhning [1999] D. Böhning,  _Computer-Assisted Analysis of Mixtures and Applications: Meta-analysis, Disease Mapping and Others_ , Monographs on Statistics and Applied Probability (Chapman & Hall/CRC, 1999).
  * McLachlan and Peel [2000] G. McLachlan and D. Peel,  _Finite Mixture Models_ , Wiley Series in Probability and Statistics (Wiley, 2000).
  * Frühwirth-Schnatter [2006] S. Frühwirth-Schnatter,  _Finite Mixture and Markov Switching Models_ , Springer Series in Statistics (Springer New York, 2006).
  * Mengersen, Robert, and Titterington [2011] K. Mengersen, C. Robert, and M. Titterington,  _Mixtures: Estimation and Applications_ , Wiley Series in Probability and Statistics (Wiley, 2011).
  * McNicholas [2016] P. McNicholas,  _Mixture Model-Based Classification_ (CRC Press, 2016).
  * McLachlan, Lee, and Rathnayake [2019] G. J. McLachlan, S. X. Lee, and S. I. Rathnayake, “Finite Mixture Models,” Annual Review of Statistics and Its Application 6, 355–378 (2019).
  * Pedregosa  _et al._ [2011] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, “Scikit-learn: Machine Learning in Python,” Journal of Machine Learning Research 12, 2825–2830 (2011).
  * Champac and Garcia Gervacio [2018] V. Champac and J. Garcia Gervacio,  _Timing Performance of Nanometer Digital Circuits Under Process Variations_, 1st ed., Frontiers in Electronic Testing No. 39 (Springer International Publishing, Cham, 2018).
  * Gregor [1969] J. Gregor, “An Algorithm for the Decomposition of a Distribution into Gaussian Components,” Biometrics 25, 79–93 (1969).
  * Xu and Mareček [2018] J. Xu and J. Mareček, “Parameter Estimation in Gaussian Mixture Models with Malicious Noise, without Balanced Mixing Coefficients,” in  _56th Annual Allerton Conference on Communication, Control, and Computing (Allerton)_ (2018).
  * Boyd and Vandenberghe [2004] S. Boyd and L. Vandenberghe,  _Convex Optimization_ (Cambridge University Press, Cambridge, UK, 2004).
  * Pearson and Henrici [1894] K. Pearson and O. M. F. E. Henrici, “Contributions to the mathematical theory of evolution,” Philosophical Transactions of the Royal Society of London. (A.) 185, 71–110 (1894).