Sarah Wassermann
email: sarah.wassermann@inria.fr

Pedro Casas
email: pedro.casas@ait.ac.at

BIGMOMAL -Big Data Analytics for Mobile Malware Detection

Keywords: • Computing methodologies → Machine learning, • Security and privacy → Network security, Mobile-Malware Detection, Big-Data Analytics, Machine Learning, High-Dimensional Data. ACM Reference format:

Mobile malware is on the rise. Indeed, due to their popularity, smartphones represent an attractive target for cybercriminals, especially because of private user data, as these devices incorporate a lot of sensitive information about users, even more than a personal computer. As a matter of fact, besides personal information such as documents, accounts, passwords, and contacts, smartphone sensors centralise other sensitive data including user location and physical activities. In this paper, we study the problem of malware detection in smartphones, relying on supervised-machine-learning models and big-data analytics frameworks. Using the SherLock dataset, a large, publicly available dataset for smartphone-data analysis, we train and benchmark tree-based models to identify running applications and to detect malware activity. We verify their accuracy, and initial results suggest that decision trees are capable of identifying running apps and malware activity with high accuracy.

INTRODUCTION

Today, more and more tasks are performed on smartphones, and they see an impressive popularity: while 472 million devices were sold in 2011 [START_REF] Goasduff | Gartner says worldwide smartphone sales soared in fourth quarter of 2011 with 47 percent growth[END_REF], this number has been more than doubled in 2014 [START_REF] Goasduff | Gartner says smartphone sales surpassed one billion units in 2014[END_REF]. Moreover, users spend most of their time browsing on their mobile phones, more than on any other device. In the United States and the United Kingdom, for instance, smartphone users spend 2-3x more hours (87 hours per month in the US and 66 hours in the UK) on their mobile device than on desktop machines (34 hours and 29 hours, respectively) [START_REF]The Communications Market Report[END_REF]. The impressive growth of mobile devices and usage has led to unprecedented increase in cellular traffic. Unfortunately, these devices increasingly become the target of malicious attacks that aim at stealing the users' private information and using it for bad-natured activities. Smart devices are an appealing target not only because of their popularity, but also because they incorporate a very large amount of sensitive information about the users, even more than a personal computer [START_REF] Suarez-Tangil | Evolution, Detection and Analysis of Malware for Smart Devices[END_REF]. For instance, the included sensors can track the users' current location and physical activities, and the embedded cameras can take pictures and record videos, even without the users' awareness.

The first piece of malware was detected before 2004 [START_REF] Suarez-Tangil | Evolution, Detection and Analysis of Malware for Smart Devices[END_REF] and the number of malware detections and their diversity grew with the expansion of the smart-device market. In particular, popular mobile operating systems such as Symbian, iOS, and Android accelerated the development of mobile malware: while only about 400 malware infections were reported between 2004 and 2007 [START_REF] Suarez-Tangil | Evolution, Detection and Analysis of Malware for Smart Devices[END_REF], this number increased by 3,325% in the last seven months of 2011 for Android alone [START_REF] Juniper | 2011 mobile threats report[END_REF]. In 2016, Kaspersky Lab found approximately 8,500,000 malicious installation packages, of which 130,000 mobile banking trojans and 260,000 mobile ransomware trojans [START_REF]Mobile Malware Evolution 2016[END_REF], and noted that a large part of malicious software is distributed via the official app stores. The year 2016 also underlined that developers of malicious apps make use of the current trends: an app called Guide for Pokémon Go including a trojan -Ztorg -made its appearance on Google Play and was installed more than 500,000 times [START_REF] Unuchek | Rooting Pokémons in Google Play Store[END_REF]. Ztorg belongs to the group of rooting malware, which gains super-user rights and is thus able to do almost everything on the infected device [START_REF]Mobile Malware Evolution 2017[END_REF]. In addition to that, the user cannot remove this kind of malware by simply performing a factory reset of her smartphone, which makes Ztorg an even more frightening threat. Even though the number of rooting-malware infections decreased in 2017, these attacks are still very popular. According to the same report, 2017 was the year of the rise of WAP and mobile banking trojans. WAP trojans steal money from the user by visiting (often without her knowing) websites with WAP subscriptions. Additionally, mobile banking trojans are now going beyond traditional financial apps and target applications designed for booking hotels, taxis, etc.

Given the relevance and sensitivity of this topic, we study the problem of malware and running-application detection in Android smartphones relying on supervised-machine-learning models. Using a large, publicly available dataset for smartphone data analysis (the SherLock data collection 1 ), we train and benchmark different supervised-machine-learning models to detect malware activity. We evaluate three different concepts, including (i) overall model performance, (ii) generalisation of the learned models across different users, and (iii) detection accuracy drift along time. Initial results suggest that tree-based models are capable of identifying malware activity with high accuracy and extremely low computational time.

The remainder of this paper is organised as follows: Section 2 discusses related work. In Section 3, we present and describe the dataset used in our study. In Section 4, we take a first approach to tackle the problem, by investigating the possibility of automatically recognising the apps running on a given device from low-level footprints, and apply the learnings and approaches later in Section 5 to automatically detect malware. Finally, Section 6 concludes this work.

RELATED WORK

Detecting mobile malware is a very challenging task and subject to significant research efforts. Zhou et al. found that the best antivirus software merely detects 79.6% of malware [START_REF] Zhou | Dissecting Android Malware: Characterization and Evolution[END_REF], which is a rather disappointing result. It is worth noting that the wide range of touch commands (such as swiping and tapping on the display) adds another layer of complexity to the process of malware detection, as this makes the exploration of all execution paths almost impossible [START_REF] Tam | The Evolution of Android Malware and Android Analysis Techniques[END_REF]. Nevertheless, the authors of Malton [START_REF] Xue | Malton: Towards on-device noninvasive mobile malware analysis for ART[END_REF] developed a system which promises to be a non-invasive Android-malware-detection tool outperforming multiple other proposals.

According to [START_REF] Suarez-Tangil | Evolution, Detection and Analysis of Malware for Smart Devices[END_REF], Android OS relies more on platform protection than on market protection. In other words, Android users can download applications even from unofficial markets, but the operating system attempts to protect the users from malicious behaviour by, for instance, applying a permission-based system (i.e., an application can only perform the actions for which the permission has been granted) or letting applications run in a sandbox (i.e., an environment isolating the app). Non-negligible efforts have been spent to detect the overprivilege problem, i.e., an application requests more permissions than it needs to work properly. Projects such as Stowaway [START_REF] Felt | Android Permissions Demystified[END_REF], its successor PScout [START_REF] Au | PScout: Analyzing the Android Permission Specification[END_REF], and DroidRisk [START_REF] Wang | Quantitative Security Risk Assessment of Android Permissions and Applications[END_REF] try to identify applications incorporating this kind of security threat. The work presented in [START_REF] Zhou | Get Off of My Market: Detecting Malicious Apps in Official and Alternative Android Markets[END_REF] uses permission-based heuristics to determine whether an application is malicious or not.

In [START_REF] Tam | The Evolution of Android Malware and Android Analysis Techniques[END_REF][START_REF] Arshad | Android Malware Detection & Protection: A Survey[END_REF], Arshad et al. and Tam et al. provide an extensive survey about malware detection and protection on Android. They present two approaches for detecting malware on Android phones: static and dynamic ones. While the static techniques are mainly based on source code analysis, dynamic methods analyse applications during their execution. For example, static techniques involve signature and permission analyses, whereas dynamic ones use anomaly detection approaches and emulation. The aforementioned works [START_REF] Felt | Android Permissions Demystified[END_REF][START_REF] Au | PScout: Analyzing the Android Permission Specification[END_REF][START_REF] Wang | Quantitative Security Risk Assessment of Android Permissions and Applications[END_REF][START_REF] Zhou | Get Off of My Market: Detecting Malicious Apps in Official and Alternative Android Markets[END_REF] fall into the category of static-analysis systems. TaintDroid [START_REF] Enck | Taint-Droid: An Information-flow Tracking System for Realtime Privacy Monitoring on Smartphones[END_REF] and Maline [START_REF] Dimjašević | Evaluation of Android Malware Detection Based on System Calls[END_REF] are two examples of systems making use of dynamic techniques. TaintDroid relies on captured network data to analyse Android applications, while Maline is a 1 http://bigdata.ise.bgu.ac.il/sherlock/ malware-detection tool which is based on machine-learning techniques to classify system calls. Another tool relying on a dynamic approach has been proposed by Afonso et al. [START_REF] Afonso | Identifying Android malware using dynamically obtained features[END_REF]: it investigates the system calls and API uses of an application to feed extracted features to a machine-learning algorithm.

To help the cybersecurity research community, several researchers make their mobile malware datasets publicly available. A good example is the SherLock dataset [START_REF] Mirsky | Sherlock vs Moriarty: A Smartphone Dataset for Cybersecurity Research[END_REF], a massive time-series dataset containing more than 600 billion samples. More precisely, they analysed the activities of Android applications in a very finegrained fashion and simulated malicious behaviours on the monitored smartphones. We use this dataset in our work and describe it in more details in Section 3. Other examples are the Device Analyzer dataset [START_REF] Wagner | Device analyzer: Large-scale mobile data collection[END_REF], a very large dataset including data instances collected from nearly 900 different devices showing the behavioural diversity among users, the LiveLab Project [START_REF] Shepard | Livelab: Measuring wireless networks and smartphone users in the field[END_REF], a dataset investigating the usage of iPhone 3GS smartphones, and the Android Malware Genome Project [START_REF] Zhou | Dissecting Android Malware: Characterization and Evolution[END_REF], a dataset encompassing more than 1,200 different Android malware samples from nearly 50 distinct malware families.

THE SHERLOCK DATASET

In the context of this work, we rely on the SherLock dataset, published by the BGU Cyber Security Research Center [START_REF] Mirsky | Sherlock vs Moriarty: A Smartphone Dataset for Cybersecurity Research[END_REF]. The SherLock dataset has been collected from January 2015 until December 2017 and includes very valuable information about the usage of smartphones. The collection was done during a long-term field trial on 50 smartphones used as primary device for different participants. More precisely, Mirsky et al. handed a Samsung Galaxy S5 to 50 clients. On these smartphones, the authors installed two applications, namely SherLock and Moriarty: while SherLock is responsible for collecting data about the smartphone (SMS, call logs, etc.), Moriarty periodically performs malicious activities, simulating the behaviour of malware found in the wild. SherLock monitors two kinds of sensors: PULL sensors which are analysed regularly (to collect, for instance, information about the installed and currently running applications, and the values of accelerometers) and PUSH sensors which are triggered when an event occurs (for example, receiving / issuing a call or a text message). The malicious actions of Moriarty are diverse and encompass activities such as contact, photo and SMS theft, phishing, and ransomware: its payload is updated every few weeks. Data collection is done at very high sample rate (as low as 5 seconds). The complete labelled dataset contains more than 600 billion data records, with a total of more than 4 TB of data. It is worth noting that the SherLock team cares a lot about the users' privacy: network identifiers such as cell tower IDs, SSIDs, and MAC addresses have been hashed; geographical locations have been anonymised through clustering techniques.

In this paper, our evaluation focuses on the data collected during the second quarter of 2016. From the complete dataset, we keep two specific feature categories: all those features related to the network traffic generated by the apps, and all those features corresponding to the footprint of the app on the CPU and internal running processes (e.g., statistics on CPUs, memory usage, Linux-level process information). The rationale is that some malware activity would be more visible at the network-traffic level, whereas some others would be better identified at the local-process level. Given the size of the data we study, we rely on big-data platforms to conduct the study. In particular, we use the Big-DAMA framework [START_REF] Casas | Network security and anomaly detection with Big-DAMA, a big data analytics framework[END_REF] to process the dataset, build machine-learning models, and evaluate their performance. In a nutshell, Big-DAMA uses Apache Spark Streaming for stream-based analysis, Spark for batch analysis, and Apache Cassandra for query and storage. Within the Big-DAMA platform, we have conceived different algorithms for network security and anomaly detection using supervised-and unsupervised-machine-learning models.

We start by doing a brief characterisation of the considered dataset snippet. During these three months, we found more than 600,000,000 records from 48 different users. Figure 1 depicts the number of different Android applications launched by each user. We can easily see that not all the users use their smartphone in the same way: while nearly 50% launched more than 200 apps, there are also clients who used less than 10 different applications. For the latter, we suppose that they did not use the Samsung smartphone on a daily basis, but only periodically for the SherLock campaign. Next, we check which applications were the most popular among our users. We consider an app as popular if the dataset contains a large number of data instances for it. Table 1 shows the ten most popular apps for Q2 2016. We note that SherLock was the app yielding by far the most samples, which comes as no surprise because it is the application responsible for the information gathering. Through this table, we can also point out that a non-negligible part of the popular applications are Android services which are automatically launched in the background by the Android OS and not by the users. For instance, SmartcardService ensures that applications can communicate with the SIM card and Google Play services allows the installed applications to access the newest features published by Google. In Table 2, we summarise the ten most popular apps excluding internal Android services and the SherLock application. According to our findings, WhatsApp and Hangouts seem to have been the favourite social applications of our users from April to June 2016.

To better understand the behaviour of Moriarty as compared to benign applications, we compare WhatsApp -the most popular social app in the dataset -against Moriarty, the malware-emulation app. In Figure 2 particular, we focus on the average CPU utilisation per hour, the average number of used threads per hour, and the sent / received traffic on a hourly basis. We observe that Moriarty uses much more CPU resources than WhatsApp: while WhatsApp used most of the time only 5% of the CPU on the monitored smartphones, Moriarty mostly needed at least 20%, which is a significant difference. However, WhatsApp seems to need more threads to run properly compared to Moriarty. When it comes to network traffic, Moriarty sent only a couple of kilobytes per hour, whereas WhatsApp sent multiple megabytes. We have a similar behaviour for the received traffic. While this is expected, as the users of WhatsApp constantly exchange messages and multimedia files, this simple comparison makes the point in the statistical differences in some of the aforementioned characteristics, which shall form the basis for the input features to the machine-learning models detecting malware activity.

IDENTIFYING RUNNING APPLICATIONS

Before setting the focus on the malware detection, we investigate the possibility of automatically identifying the different running Android applications from the set of input features available in the dataset. The goal is to blindly differentiate among different running apps (i.e., without gathering the process name), so as to better understand if it would be possible to later on identify malware activity using the same approach.

We consider as input a set of 45 SherLock-monitored smartphone features, reflecting the behaviour of the different apps, and build supervised-machine-learning models to automatically classify them. We particularly focus on the usage of decision trees, based on their a-priory excellent performance in similar network-measurement analysis problems [START_REF] Wassermann | NETPerfTrace: Predicting Internet Path Dynamics and Performance with Machine Learning[END_REF]. We proceed as follows: for each user, we randomly choose 20% of the application samples as training set and the remaining 80% for testing. We obtain almost perfect results for the entire set of Q2 users: for all of them, our model correctly identified the running applications in about 99% of the cases.

Feature Selection

Next, we investigate which features are the most discriminative ones when it comes to the identification of applications. We rely on an embedded method, based on decision trees and the notion of variable importance [START_REF] Geurts | Supervised learning with decision treebased methods in computational and systems biology[END_REF]. With this kind of method, the most relevant features are determined while building the prediction model. The ten most important features are reported in Table 3. We note that the most discriminative features are related to the CPU, threads, and memory. Moreover, the information related to the importance of the Android application -i.e., whether it is running as a service, in the background, or in the foreground -helps for the identification. Overall, our results show that low-level information is very useful for this kind of prediction task. As an additional analysis, we evaluate the application-identification power of a decision tree for one user when relying only on the five and ten most important features. For this, we split the dataset as previously, i.e. 20% of the data is considered as the training set and the remaining 80% as testing set. Our results reveal that the decision-tree model still works extremely well even when using only 5 or 10 features instead of 45: in both cases, the identification was correct for about 99% of the testing records. It thus seems that only a small amount of data is required to obtain highly accurate prediction performance, which has a paramount advantage in case the app identification needs to be done in an online context. Indeed, when working with a restricted set of features, classifiers can work more efficiently and provide an answer quickly. Furthermore, gathering only little data implies that we do not have to overload the smartphone for data collection -only a few system calls are necessary.

Temporal App Identification

Finally, we are interested in the quality of the application identification on a temporal basis. More precisely, we are curious about the accuracy of the prediction model in case the training phase takes place in a different time span. In particular, for each user, we train a decision tree on data collected during the first week and try to identify the applications run during the following weeks. Figure 3, depicting the results for one user, reveals that the overall identification accuracy decreases with respect to our previous evaluation, suggesting that restricting the training set to the samples collected during the first week does not capture as well the users' smartphone habits as when the training set is randomly sampled. Furthermore, the plot shows that the precision of our model decreases over time: while almost 80% of the applications were correctly identified in the first prediction week, this was the case for only about 62% in the last one. The outcome of this analysis for the other users is very similar.

MALWARE DETECTION

Our next goal is to detect whether an application is malicious or not, still working with the data collected during the second quarter of 2016. As presented in Section 3, Moriarty simulates bad-natured behaviour and is thus the target of the detection model. We are now trying to answer the following question: "Is the running application Moriarty or not?", which boils down to a binary classification problem. We use the same feature set and prediction model as for the application-identification task. However, as the execution of the Moriarty application can be considered as a rare event (the dataset contains millions of non-Moriarty samples versus only a few thousands of Moriarty records), we randomly choose 40% of samples for the testing set compared to 20% in the context of the application identification. Figure 4 depicts the obtained recall for the Moriarty detection for each user. The recall is defined as #true positives #true positives + #f alse neдatives , and a high value thus reveals that the number of missed malware instances (the false negatives) is low. This is a very important property in our context, as false negatives (i.e., considering an application as benign although it is malicious) are a considerable security threat for the user. Moreover, we achieve an extremely low false positive rate (below 1% for each user), indicating that we almost never classify a benign application as malicious. For the sake of illustration, Table 4 shows the confusion matrix obtained for one of our users.

Feature selection

Similar to the identification of applications, we are interested in the most important features for detecting malicious apps. The ten most relevant features (determined using the same technique as in Section 4) are summarised in Table 3. We note that, compared to the application-identification task, a lot of these features are related to the memory (five out of ten). Once again, low-level metrics are very useful indicators. As a final step, and similarly to Section 4, we evaluate our decision-tree model for this detection problem with the top 5 and top 10 features. The results are excellent even with these reduced feature sets. Indeed, in both cases, the recall is over 99%.

Temporal Malware Detection

Next, we carry out the temporal analysis for the Moriarty detection. We proceed in the same way as for our previous task, i.e., we train for each user the detection model on data gathered during the first week of the second quarter and test it in the remaining weeks.

Given the major imbalance in classes between benign and Moriarty instances, we resort to standard bootstrapping to rebalance classes for learning purposes. In particular, we randomly oversample the Moriarty class. Figure 5 depicts the recall obtained per week for one user. As we can easily observe, the results are very poor and detection performance rapidly degrades after the 3rd week. Two issues could be the cause of this performance: firstly, the extremely low number of gathered Moriarty samples in the training week, and secondly, a concept drift; as stated in Section 3, the Moriarty application periodically changes its behaviour, which makes its detection very challenging when the training set does not include samples of all the patterns. Indeed, during the second quarter of 2016, Moriarty changed twice its malicious behaviour: it started with SMS / banktheft malware activity, moved on to phishing, and finished with simulating adware [START_REF] Mirsky | Sherlock vs Moriarty: A Smartphone Dataset for Cybersecurity Research[END_REF]. This prevents the decision tree from accurately grasping the behaviour of the malicious Android application. Again, the results for the other users are highly similar. We tested several other models for this specific task, without any notable detection improvements.

Malware Detection across Multiple Users

Finally, we evaluate our detection model on data collected from multiple users. More precisely, we train our decision tree on data gathered from three Q2 users and test it on data belonging to three other ones. The resulting confusion matrix is shown in Table 5. Results are again quite disappointing, suggesting that the particularities of the ways users use their smartphone play a more important role than the set of features used as input by the model. Indeed, each and every user uses her smartphone differently and it is thus very challenging to build a model from multiple users that can be accurately applied to other ones. In other words, malware detection based on behavioural features from a community, i.e., features describing how the smartphone is used, is a very difficult task.

Predicted positive Predicted negative

CONCLUDING REMARKS AND PERSPECTIVES

In this paper, we tackled two prediction tasks in the domain of smartphone security by relying on the SherLock dataset: the identification of running Android-applications and the detection of malware. For both of these classification targets, we obtained highly promising results in case the training sets are generated randomly. Nevertheless, the Android application identification and malware detection on a temporal basis as well as when relying on data instances from multiple users remains a challenging task. There are two main takeaways from this study: firstly, our results suggest that malware-detection models using the set of features provided by the SherLock dataset can operate with very high performance only in those cases where models are tailored on a per-user basis, taking into account the statistical behaviour of each user. Secondly, it seems hard to build a model that generalises well and can cover multiple users and longer time spans, additionally due to concept drift occurrences. We are therefore currently investigating stream-based machine-learning approaches for the malware-detection task, in which we periodically retrain the underlying model when a concept drift and / or a detection performance degradation is detected.

Figure 1 :

 1 Figure 1: Number of different applications launched during Q2 2016 by each user.

  Avg. CPU utilisation / hour for WhatsApp. (b) Avg. # used threads / hour for WhatsApp. (c) # sent megabytes / hour by WhatsApp. (d) # received megabytes / hour by WhatsApp. Avg. CPU utilisation / hour for Moriarty. (f) Avg. # used threads / hour for Moriarty. (g) # sent kilobytes / hour by Moriarty. (h) # received megabytes / hour by Moriarty.

Figure 2 :

 2 Figure 2: Characterisation of the WhatsApp and Moriarty Android applications. Each curve on the CDFs corresponds to one user and each point of the curve to one hour.

Figure 3 :

 3 Figure 3: App identification accuracy when training and testing on a weekly basis.

Figure 4 :

 4 Figure 4: Recall obtained while detecting Moriarty for each Q2 user.

Figure 5 :

 5 Figure 5: Recall obtained while detecting Moriarty on a weekly basis for one Q2 user.

Table 1 : Top 10 popular apps in Q2 2016.

 1 , we plot several characteristics of both of them. In

	Rank	Application name
	#1		SherLock
	#2		Google Play Services
	#3		Chrome
	#4		SmartcardService
	#5		Unified Deamon
	#6		SecurityManagerService
	#7		WhatsApp
	#8		S Finder
	#9		Samsung Push Service
	#10		Context Service
	Rank		Application name
	#1		Chrome
	#2		Google App
	#3		WhatsApp
	#4		S Finder
	#5		S Health
	#6		Hangouts
	#7	Samsung text-to-speech engine
	#8		Geo News
	#9		Peel Smart Remote
	#10		Contacts

Table 2 : Top 10 popular apps in Q2 2016 (Android services and SherLock excluded).

 2

Most relevant features for identifying running Android applications and detecting malware.

	Top 10 features	Application identification	Malware detection
	#1 feature	virtual memory size	shared dirty pages used by Dalvik heap
	#2 feature	number of threads	number of minor page faults
	#3 feature	CPU utilisation	time process scheduled in user mode
	#4 feature	process priority level (foreground, background, service, sleeping, etc.)	proportional set size for Dalvik heap
	#5 feature	ordering within a particular priority category	time process-children scheduled in user mode
	#6 feature	process life time	virtual memory size
	#7 feature	time process scheduled in user mode	time process scheduled in mode
	#8 feature	time process-children scheduled in user mode	process ID
	#9 feature	number of minor page faults	process life time
	#10 feature	number of private dirty pages used by everything else than the native	number of private dirty pages used by the native heap
	3:

Table 4 : Confusion matrix obtained while detecting Mori- arty for one Q2 user. Predicted positive Predicted negative

 4

	True positive	14,275,309	155
	True negative	117	35,076
	True positive		16,784
	True negative	35,401	3,472

Table 5 : Confusion matrix obtained while detecting Mori- arty for a couple of Q2 users when training on data instances from other clients.

 5

partially funded by the Vienna Science and Technology Fund (WWTF) through project ICT15-129, "BigDAMA". * This work has been carried out while the author was an MSc. student at the University of Liège.


Lizheng Zhang
email: lizhengz@cae.wisc.edu

Weijen Chen
email: weijen@cae.wisc.edu

Yuhen Hu
email: hu@engr.wisc.edu

Chung-Ping Chen
email: chen@engr.wisc.edu

Statistical Timing Analysis with Extended Pseudo-Canonical Timing Model

State of the art statistical timing analysis (STA) tools often yield less accurate results when timing variables become correlated due to global source of variations and path reconvergence. To the best of our knowledge, no good solution is available dealing both types of correlations simultaneously.

In this paper, we present a novel extended pseudo-canonical timing model to retain and evaluate both type of correlation during statistical timing analysis with minimum computation cost. Also, an intelligent pruning method is introduced to enable trade-off runtime with accuracy.

Tested with ISCAS benchmark suites, our method shows both high accuracy and high performance. For example, on the circuit c6288, our distribution estimation error shows 15× accuracy improvement compared with previous approaches.

I. INTRODUCTION

It is well-known that the timing performance of future generations of deep-submicron micro-architecture will be dominated by several factors. IC manufacturing process parameter variations will cause device and circuit parameters to deviate from their designed value. Low supply voltage for lowpower applications will reduce noise margin, causing increased timing delay variations. Due to dense integration and non-ideal on-chip power dissipation, rising temperature of substrate may lead to hot spot, causing excessive timing variations.

Classical worst case timing analysis produces timing predictions that are often too pessimistic and grossly conservative. On the other hand, statistical timing analysis (STA) that characterizes timing delays as statistical random variables offers a better approach for more accurate and realistic timing prediction.

In literatures, there are two distinct approach for STA: path based STA and block based STA. The fundamental challenge of the path based STA [START_REF] Liou | False-path-aware statistical timing analysis and efficient path selection for delay testing and timing validation[END_REF]- [START_REF] Agarwal | Statistical delay computation considering spatial correlations[END_REF] is its requirement to select a proper subset of paths whose time constraints are statistically critical. This task has a computation complexity that grows exponentially with respect to the circuit size, and hence can not be easily scaled to handle realistic circuits.

This potential difficulty has motivated the development of block base STA [START_REF] Agarwal | Statistical timing analysis using bounds and selective enumeration[END_REF]- [START_REF] Chang | Statistical timing analysis considering spatial correlations using a single pert-like traversal[END_REF] that champions the notion of progressive computation. Specifically, statistical timing analysis is performed block by block in the forward direction in the circuit timing graph without looking back to the path history. As such, the computation complexity of block based STA will grow linearly with respect to the circuit size. To even further speed up the computation, Gaussian assumption has been widely adopted( [START_REF] Visweswariah | First-order parameterized block-based statistical timing analysis[END_REF], [START_REF] Devgan | Block-based static timing analysis with uncertainty[END_REF], [START_REF] Chang | Statistical timing analysis considering spatial correlations using a single pert-like traversal[END_REF]) with small accuracy penalty, and all internal timing random variables in a circuits are forced to follow the Gaussian distribution.

However, to realize the full benefit of block based STA, one must solve a difficult problem that timing variables in a circuit could be correlated due to either global variation ( [START_REF] Visweswariah | First-order parameterized block-based statistical timing analysis[END_REF], [START_REF] Agarwal | Statistical timing analysis for intra-die process variations with spatial correlations[END_REF], [START_REF] Chang | Statistical timing analysis considering spatial correlations using a single pert-like traversal[END_REF]) or path reconvergence( [START_REF] Agarwal | Statistical timing analysis using bounds and selective enumeration[END_REF], [START_REF] Devgan | Block-based static timing analysis with uncertainty[END_REF]). As illustrated in the left hand side of Figure 1, global correlation refers to the statistical correlation among timing variables in the circuit due to global variations such as inter-or intra-die spatial correlations, same gate type correlations, temperature or supply voltage fluctuations, etc. Path correlation, illustrated in the right hand side of Figure 1, refers to the correlation resulting from the phenomenon of path reconvergence, that is, timing variables may share a common subset of gate or interconnect along their path histories. The importance of the path correlation comes from the fact that each gate in the circuit will have some local variations which are independent to other gates in the circuit. These local variations will propagate towards the circuit output and cause additional correlations. However, the correlations caused by sharing some of these local variations because of path reconvergence, cannot be correctly captured by any algorithm that deals with global variations only.

Several preliminary solutions have been proposed to deal with these correlations. In [START_REF] Visweswariah | First-order parameterized block-based statistical timing analysis[END_REF], [START_REF] Agarwal | Statistical timing analysis for intra-die process variations with spatial correlations[END_REF], [START_REF] Chang | Statistical timing analysis considering spatial correlations using a single pert-like traversal[END_REF] the dependence on global variations is explicitly represented using a canonical timing model. However, none of these approaches has taken into account the path correlations. In [START_REF] Devgan | Block-based static timing analysis with uncertainty[END_REF], a method based on common node detection is introduced to deal with the path correlations. However, this method does not address the issue of dependence on global variations.

In this paper, we present a systematic STA solution that takes into account correlations caused by both global variations and path reconvergence. Specifically,

• We extend the commonly used canonical timing model to represent all timing variables in the circuit as a weighted linear combination of a set of Gaussian random variables. A variation vector, consisting of all the weights, is then used to explicitly represent both global and path correlation information. • We develop a pseudo-canonical timing model which relax the requirement of the independence between the global variation sources involved in the original canonical timing model. With such, the expensive principle component analysis is avoided and total computation time is saved. • We further explore the sparse structure of the variation vector and develop a flexible vector format so that the non-significant entries of the variation vector are dynamically dropped during computation. According to simulations on ISCAS circuits, this technique significantly curtails the amount of storage and computation required for our method. Since min(X, Y ) = -max(-X, -Y ), in the interests of brevity, in the rest of this paper, we will only discuss the MAX operator, with the understanding that the same results can be easily adapted to the MIN operator.

The rest of the paper is organized as following: In section II, previous block based STA methods are reviewed briefly; Section III describes the vectorized timing format and a theorem used for correlation decompose; Section IV is the detailed algorithm and technique to reduce computation complexity. Section V presents a real implementation of our method in C/C++ and the testing result with ISCAS85 benchmark suites; Section VI gives the conclusions.

II. CANONICAL TIMING MODEL

In timing analysis field, the circuit is modeled as timing graph, where nodes are used to represent the gate/wire in the circuit and there will be some delay, called node delay, associated with them. Signals propagate through these nodes will add their delays into its arrival time. Node delay and arrival time are usually called timing variables of the circuit.

Different from classical timing analysis, the statistical timing analysis models timing variables as random variables, which are characterized by its probability density function(p.d.f.) or cumulative distribution function(c.d.f.). The purpose of statistical timing analysis is then to estimate the arrival time distribution at the primary output of the circuits knowing input arrival time distributions and all internal node delay distributions. This is accomplished through two operators [START_REF] Agarwal | Statistical timing analysis using bounds and selective enumeration[END_REF]:

• ADD: When an input arrival time X propagates through a node delay Y , the output arrival time will be Z = X+Y • M AX: When two arrival times X and Y merge in a node, a new arrival time Z = max(X, Y ) will be computed before the node delay is added. [START_REF] Visweswariah | First-order parameterized block-based statistical timing analysis[END_REF], [START_REF] Agarwal | Statistical timing analysis for intra-die process variations with spatial correlations[END_REF], [START_REF] Chang | Statistical timing analysis considering spatial correlations using a single pert-like traversal[END_REF] proposed a canonical delay model to address the node delay correlations through sharing global variations. In particular, they model each of the node delay as a summation of three terms:

n i = µ i + α i R i + j β i,j G j ( 1
)
where n i (i = 1, 2, ...) are random variables corresponding to the the i th node delay in the timing graph; µ i is the expected value of n i ; R i , (named node variation), is a zero-mean, unity variance Gaussian random variable representing the localized statistical uncertainties of n i ; G j represents the j th global variation, and is also modeled as a zero-mean, unity variance Gaussian random variable; {R i } and {G j } are additionally assumed to be mutually independent; the weight parameters α i (named node sensitivity) and β i,j (named global sensitivities) are deterministic constants, explicitly expressing the amount of dependence of n i on each of the corresponding independent random variables. With this canonical representation, the correlation (covariance) between any two node delays, n i and n k , can be easily evaluated.

cov(n i , n k ) = E{(n i -µ i )(n k -µ k )} = j β i,j β k,j (2) Note that random variables {R i , R k , G j (j = 1, 2, ...)} are mutually independent.
This canonical timing model is very good for node delay(gate delay or interconnect delay) modeling. But it is not sufficient to model the arrival time. When an arrival time is somehow mapped into this canonical timing model as done in [START_REF] Chang | Statistical timing analysis considering spatial correlations using a single pert-like traversal[END_REF], some important path information will be dropped and the accuracy will be compromised.

For example, referring back to Figure 1(b), arrival time X and Y share the common history of node p. So between them, there will be some extend of correlation caused by the node p's local variation R p . But X and Y have only one term of local variation of R X and R Y respectively in their canonical representations. And according to the assumption from the canonical timing model, R X and R Y are assumed to be independent to each other, so the correlation caused by R p in node p is incorrectively ignored.

Another significant disadvantage of the original format of the canonical model is that it requires the independence between the global variation sources G i . This requirement, as stated in [START_REF] Chang | Statistical timing analysis considering spatial correlations using a single pert-like traversal[END_REF], will be satisfied by Principle Component Analysis(PCA) before the timing analysis. But we think it is not reasonable to exclude the computation complexity of PCA from the overall performance evaluation of the statistical timing analysis because the complexity of the PCA grows as quickly as O[M 3 ] where M is the total number of global variation sources considered. Due to the important spatial correlation between global variation sources, the number of global variation sources have to be included in the PCA will easily go up to thousands and so that PCA itself will be a very timing consuming step.

For example, to model the spatial correlation, the chip is usually partitioned into grids and each grid will be associated with m global variations. So the total number of global variation sources will be proportional to the number of grids used to represent the chip area. If the chip is partitioned into 10×10 grid, the total number of global variations will be M = 100m. This number maybe tolerable but the coarse grid may not be sufficient to accurately model the spatial correlation of global variations between different grid locations. If we partition the chip with 100 × 100 to get more accurate spatial correlation model, the total number of global variations will now be M = 10000m which is beyond the reasonable range for cases where a fast PCA is demanded. To be a fact, for a real chip with dimension of 10mm × 10mm, 100 × 100 grids will only have grid cell with size of 100µm × 100µm which is not yet small enough for high quality spatial correlation modeling.

III. EXTENDED PSEUDO-CANONICAL TIMING MODEL

The canonical timing model [START_REF] Visweswariah | First-order parameterized block-based statistical timing analysis[END_REF], [START_REF] Agarwal | Statistical timing analysis for intra-die process variations with spatial correlations[END_REF], [START_REF] Chang | Statistical timing analysis considering spatial correlations using a single pert-like traversal[END_REF] is a powerful tool to represent the numerous timing variables for a given circuit. However, as pointed out in the previous section, in its original format, it can only handle node delay correlations caused by global variations. Furthermore, the original format of canonical timing model requires a computationally expensive PCA which may prohibit it from using in cases where large number of global variations are correlated with each other. In this work, we propose an extended pseudo-canonical timing model(EPCT) that is capable of capture all the correlation between any pair of timing variables in the circuit be it a node delay or an arrival time. And also it removes the requirement of a PCA procedure for its validation.

A. Extended Pseudo-Canonical Timing Model

We notice an important fact that it is not necessary to have the canonical format to be able to evaluate the correlations caused by global variation sources. So we propose a Pseudo-Canonical Timing Model for node delay as following:

n i = µ i + α i R i + j β i,j G j (3)
which, compared with the commonly used canonical model in equation ( 1), relaxes the requirement of the independence of the global variations and so that the global variations G j may be CORRELATED.

With such pseudo-canonical timing model, the correlation between two node delays can be evaluated as:

cov(n 1 , n 2 ) = i j β 1,i β 2,j cov(G i , G j ) (4)
Assume that there are N nodes and M global variations in the timing graph, if every node delay can be modeled by the pseudo-canonical format of Equation ( 3), then every timing variable, including all the node delays and arrival times will then have a extended pseudo-canonical timing(EPCT) model as:

X = µ X + N i=1 α X,i R i + M j=1 β X,j G j ( 5
)
where R i and G j are independent to each other although there will possibly correlations between G j s. With this equation, both global and path correlations can be handled elegantly. More specifically, global variations are represented by the set of global sensitivity terms {β X,j }, and dependence on path history are represented by non-zero node sensitivity terms α X,k .

B. Variation Vector

Equation ( 5) can be rewritten in a compacted vector format as

X ∼ L(µ X , α X , β X ) = µ X + α * X r + β * X g (6)
where "*" means transpose and

r ≡ [R 1 , • • • , R N ] * ∼ N (0, I) g ≡ [G 1 , • • • , G M ] * ∼ N (0, Σ g ) (7)
are mutually independent local variation vector and global variation vector respectively. 0 is a zero vector and Σ g = E{gg * } is the covariance matrix of global variations and generally not equal to the unit matrix I due to the potential correlations existing between global variations. α

X = [α X,1 , α X,2 , ..., α X,N ] * and β X = [β X,1 , β X,2 , ..., β X,M ] * are deterministic vectors called Variation Vector(v.v.) of X.
Authors in [START_REF] Chang | Statistical timing analysis considering spatial correlations using a single pert-like traversal[END_REF] proves the correlation evaluation formula between timing variables represented by the canonical timing model of equation [START_REF] Liou | False-path-aware statistical timing analysis and efficient path selection for delay testing and timing validation[END_REF]. We here prove a similar formula for correlation evaluation between time variables expressed with the EPCT model as equation ( 5) or [START_REF] Visweswariah | First-order parameterized block-based statistical timing analysis[END_REF].

Theorem 1:

Given timing variables X ∼ L(µ X , α X , β X ) and Y ∼ L(µ Y , α Y , β Y )
, the correlation between them can be evaluated as:

cov(X, Y ) = α * X α Y + β * X Σ g β Y (8) Proof: By definition: cov(X, Y ) = E{(X -µ X )(Y -µ Y )} = cov(α * X r, α * Y r) + cov(α * X r, β * Y g) +cov(α * Y r, β * X g) + cov(β * X g, β * Y g) = E{α * X rr * α X } + E{β * X gg * β X } = α * X α Y + β * X Σ g β Y
where the independence between r and g is applied.

For the variance of a time variable, it is easy to get: Corollary 1: Given timing variable X ∼ L(µ X , α X , β X ), its variance is:

σ 2 X = α * X α X + β * X Σ g β X ( 9
) This corollary is actually the special case when X = Y of theorem 8.

IV. PROPAGATING MEAN AND VARIATION VECTOR

In a timing graph, the mean and variation vectors of a node delay is obtained from technology extraction. A STA algorithm, instead, will take those node's means and variation vectors as its input and calculate the mean and variation vector for all arrival times in the entire circuit.

A. Exploration of Sparsity

Since there are N nodes and M global variations in a timing graph, the length of the local variation α will be N and the length of global variation vector β will be M which is the same as the rank of the correlation matrix Σ g . So for every step of arrival time propagation in the timing analysis, the computation complexity will be O(N + M 2 ). And the overall timing analysis computation will have complexity of O(N (N + M 2 )) ≈ O(N 2 + N M 2 ) since there will be N steps of arrival time propagations.

Although the correlation matrix Σ g will be dense and have rank of M which is the same as the total number of global variation sources. However, while working with benchmark circuits, we noticed that many components in the variation vectors have very small values, indicating that their contributions to the overall correlation evaluation is insignificant. By setting these small coefficient to zero, the variation vector will become a sparse vector that contains many zero components.

Motivated by this observation, we developed a novel technique called the flexible vector format to exploit the sparsity of the variation vector. For this purpose, a drop threshold is selected so that if α X,i or β X,j is smaller than this threshold, it is deemed to have small value and will be placed into a drop candidate pool to be pruned from the variation vector representation.

However, dropping α X,i or β X,j with small magnitude is the same as applying truncation to the variation vector. In subsequent computations, the quantization error may accumulate, causing non-negligible error. This is a problem that can not be overlooked for large circuits. Our solution to this problem is to lump those components in the drop candidate pool into a single correction term

x pool = x 2
dropped Components (10)

B. Complexity and Path Correlation Length

Using this drop and pool mechanism, the actual number of non-zero terms in β will be much smaller than M . Assuming there are m global variation sources for a single node, the average number of non-zero terms in β will be either (1) m if the time variable is a node delay; or (2) M C = m × q << M if the time variable is an arrival time and q is proportional to the logic depth of the circuit. So, using equation ( 8) to evaluate the correlation between time variables, the complexity for a single evaluation will have complexity of O[M 2 C + Γ] where the average number of non-zero terms in local variation vector α is Γ.

Using this drop and pool mechanism, what is really dropped in the local variation vector α during computation is then the path correlation. So the length of the local variation vector actually gives a good indication to the extent of path correlation in the circuit. The path correlation length(Γ) of the circuit is then defined as the average length of the pruned local variation vectors for a given drop threshold.

Theoretically speaking, the path correlation length will be the number of critical nodes which are in the statistically critical paths. Specifically, If a node is not in any statistically critical paths, its variation will be automatically dropped. On the other hand, if the node is in one critical path but is not statistically important, it will be dropped too. So by studying the variation vector in the circuit output, it is easy to know which path is critical and which node is critical for circuit performance and this is very important information to help designer improve the circuit design.

In real circuits, usually only a few paths are statistically critical and so that only a few nodes in the circuit will survive the variation vector propagation. So the computation complexity of our method will be O

[(Γ + M 2 C ) • N ].
Since Γ << N and M C << M , the total complexity will then be O[N ] and a significant reduction of computation and storage is achieved with the drop and pool mechanism. Compared with the PCA-based approach, the complexity saving is even larger since there the complexity will be 3 + N ] even the drop and pool mechanism is used. So the extended pseudo-canonical timing model actually results in a big saving in computation complexity when the total number of global variations is large.

O[M 3 + (M C + Γ)N ] ≈ O[M

V. SIMULATION RESULTS AND DISCUSSIONS

The above described algorithm has already been implemented in C/C++ and tested by ISCAS85 benchmark circuits.

Before testing, however, all benchmark circuits are remapped into a library which has gates of not, nand2, nand3, nor2, nor3 and xor/xnor. All library gates are implemented in 0.18µm technology and their delays are characterized by Monte Carlo simulation with Cadence tools assuming all variation sources, either process variations or operational variations, follow Gaussian distribution.

For illustration purpose, only three parameter variation are considered global: channel length(L), supply voltage(Vdd) and temperature(T). All other variation sources, specified in the 0.18µm technology file, are assumed to be localized in the considered gate only.

A. Accuracy and Performance

Monte Carlo simulation results with 10,000 repetitions are used as "Golden Value" for each benchmark circuit. Each repetition is a process of static timing analysis by fixing global and node variation into a set of randomly sampled values. The global variations are sampled once for each repetition while node variation for each gate is newly sampled every time when the gate is computed.

To elaborate the importance of including the path reconvergence correlation, A STA methods called CTM is implemented with canonical timing model where no path correlations are considered.

Table I summarizes the arrival time distribution parameters at the primary output of each testing circuit from Monte Carlo(M.C.), CTM and our method of EPCTM which is implemented with extended pseudo-canonical timing model. µ and σ are mean and standard variation of the distribution. τ 97 = µ + 2σ is the delay estimation at confidence level of 97%. The accuracy of STA methods compared with Monte Carlo method, is evaluated in Table II Table II shows that CTM have significantly larger error in mean estimation than EPCTM. This is reasonable because CTM will overestimate the mean at every MAX operation due to smaller correlation considered and this over estimation is accumulated through distribution propagation. It is also interesting to notice that CTM and EPCTM give similar accuracy in variance estimation. This is possibly because of the fact that the variance is dominated by global variation in the tested cases.

Of course, Monte Carlo simulation gives the best STA results but with big runtime penalty. EPCTM runs order-ofmagnitude faster but can provide both mean and variance estimation almost as accurate as Monte Carlo does if most of the path correlations are considered as the cases of EPCTM shown in Table I and II.

To further elaborate the accuracy of EPCTM, Figure 2 shows the p.d.f. and c.d.f. for circuit c6288 from three methods: Mont Carlo and two methods of EPCTM and CTM. Apparently enough, EPCTM shows excellent accuracy since most path correlation is considered.

B. Performance and Path Correlation Length

It has been mentioned in Section IV-B that path correlation length (Γ) is an interesting macro property of the simulated circuit and gives a good indication of the extent of the path correlation existing in that circuit. For the above ISCAS III where the run time improvement is also shown when EPCTM is compared with Monte Carlo.

From Table III, we can firstly conclude that the correlation length Γ is much smaller than the circuit size and basically independent on the circuit size since it remains about 10 -20 when circuit size changes dramatically. This observation helps the conclusion we made before about the complexity reduction of our method by using the technique of flexible vector format.

Secondly, the only exceptional high path correlation length among the tested circuits happens with the circuit c6288 which is known as a 16-bit array multiplier. Since there are large amount of equal delay paths in the circuit, large To study the relationship between path correlation length and the accuracy of the STA method, An experiment is conducted for circuit c6288 and results are shown in figure 4 where the error in τ 97 and path correlation length are both plotted against the drop threshold. It is clear that the path correlation length drops sharply when the drop threshold changes slightly from zero and maintain almost constant after that. But the error changes readily when drop threshold changes. This phenomenon proves the efficiency of the drop mechanism introduced in this work since it means we can sacrifice very little accuracy to gain very significant reduction in the path correlation length and so that save significant amount of CPU time since the run time of the proposed STA method is proportional to the path correlation length.

VI. CONCLUSIONS

This paper presents a novel method for block-based statistical timing analysis. Applying the generally accepted Gaussian assumption, we firstly disclose that the MAX operation can be approximated by linear supposition of its inputs. Secondly we extend the commonly used canonical timing model into a (δτ 97 ) when drop threshold changes vectorized format, variation vector. We also disclose a novel method to decompose correlated timing variables into independent ones to simplify computation. With these theoretical progress, we are able to evaluate and propagate the global and path correlation systematically in the circuit timing graph. We also design a novel algorithm which treat both global and path correlation simultaneously and systematically. This algorithm, with the help with a new flexible vector format achieves high accuracy and high performance at the same time as tested by ISCAS circuits and compared with Monte Carlo results.

Fig. 1 :

 1 Fig. 1: Global Correlations (left) and Path Correlation(right)

  Fig. 2: p.d.f. and c.d.f. comparison for c6288 from three methods

Fig. 3 :

 3 Fig. 3: Run time Complexity of EPCTM Shown in the figure 3 is the run time complexity of the proposed timing algorithm EPCTM where the run time and circuit size of all circuits except c6288 are shown. From the figure, it is clear that the run time is almost linear with respecting to the circuit size even when the circuit size changes dramatically. This result clearly demonstrates our complexity discussion in section IV-B.To study the relationship between path correlation length and the accuracy of the STA method, An experiment is conducted for circuit c6288 and results are shown in figure4where the error in τ 97 and path correlation length are both plotted against the drop threshold. It is clear that the path correlation length drops sharply when the drop threshold changes slightly from zero and maintain almost constant after that. But the error changes readily when drop threshold changes. This phenomenon proves the efficiency of the drop mechanism introduced in this work since it means we can sacrifice very little accuracy to gain very significant reduction in the path correlation length and so that save significant amount of CPU time since the run time of the proposed STA method is proportional to the path correlation length.

Fig. 4 :

 4 Fig. 4: Path correlation length(Γ) and Error in 97% delay (δτ 97 ) when drop threshold changes

  .

		STA	CPU	Delay Distribution[ps]
	Circuit	Method	Time[s]	µ	σ	τ 97
		M.C.	6.449	1288.8	219.3	1727.5
		CTM	0.010	1348.6	216.0	1780.7
	c432	EPCTM	0.030	1299.0	220.4	1739.9
		M.C.	8.182	1073.6	178.9	1431.4
		CTM	0.010	1125.4	178.3	1482.0
	c499	EPCTM	0.030	1084.8	180.5	1445.8
		M.C.	14.831	1445.4	266.3	1977.9
		CTM	0.010	1463.1	264.2	1911.5
	c880	EPCTM	0.050	1447.6	264.9	1977.3
		M.C.	19.007	1445.4	251.4	1948.3
		CTM	0.010	1529.4	249.0	2027.4
	c1355	EPCTM	0.071	1460.9	250.7	1962.3
		M.C.	35.801	1828.2	327.3	2482.8
		CTM	0.030	1881.7	326.5	2534.7
	c1908	EPCTM	0.150	1841.9	328.4	2498.6
		M.C.	72.163	2097.0	382.9	2862.8
		CTM	0.050	2161.8	379.7	2921.2
	c2670	EPCTM	0.181	2104.4	382.9	2870.1
		M.C.	84.020	2747.2	498.8	3744.8
		CTM	0.050	2850.3	500.6	3851.5
	c3540	EPCTM	0.240	2752.3	502.1	3756.5
		M.C.	140.832	2399.3	441.7	3282.6
		CTM	0.080	2474.1	441.3	3356.7
	c5315	EPCTM	0.641	2404.8	442.2	3289.2
		M.C.	114.235	6740.1	1286.8	9313.6
		CTM	0.070	7290.8	1273.1	9836.9
	c6288	EPCTM	5.198	6775.9	1275.1	9326.2
		M.C.	202.972	1911.7	348.7	2609.0
		CTM	0.110	1974.3	352.5	2679.3
	c7552	EPCTM	0.571	1916.6	353.8	2624.2

TABLE I :

 I Testing Results for ISCAS Benchmarks

  comparison for c6288 from three methods

	Circuit		Mean Error(δµ)	Variance Error(δσ)
		CTM	EPCTM	Improve	CTM	EPCTM
	c432	4.64%	0.79%	5.9x	1.50%	0.50%
	c499	4.82%	1.04%	4.6x	0.34%	0.89%
	c880	1.22%	0.15%	8.1x	0.79%	0.53%
	c1355	5.81%	1.07%	5.4x	0.95%	0.28%
	c1908	2.93%	0.75%	3.9x	0.24%	0.27%
	c2670	3.09%	0.35%	8.8x	0.84%	0.00%
	c3540	3.75%	0.19%	19.7x	0.36%	0.66%
	c5315	3.12%	0.23%	13.6x	0.09%	0.11%
	c6288	8.17%	0.53%	15.4x	1.06%	0.65%
	c7552	3.27%	0.25%	13.1x	1.09%	1.46%

TABLE II :

 II

Distribution Errors circuits, the path correlation length(Γ) at drop threshold of 1% is summarized in Table

TABLE III

 III

				: Path Correlation Length and Runtime
				Improvement over Monte Carlo
	path correlation length is natural: Few node variation can be
	dropped due to the equal importance.
					Complexity of EPCTM
		1
		0.9
		0.8
		0.7
	Run Time[s]	0.4 0.5 0.6
		0.3
		0.2
		0.1
		0	0	1000	2000	3000	4000	5000	6000
						Circuit Size

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)


Praveen Ghanta
email: ghanta@ece.arizona.edu

Sarma Vrudhula

Rajendran Panda
email: rajendran.panda@freescale.com

Janet Wang

Stochastic Power Grid Analysis Considering Process Variations

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

One of the most difficult and important challenges posed by deep sub 100 nm IC technology is the increasing uncertainty in the performance of CMOS circuits due to variations in the fabrication process [START_REF] Mehrotra | A methodology for modeling the effects of systematic within-die interconnect and device variation on circuit performance[END_REF][START_REF] Boning | Models of process variations in device and interconnect[END_REF]. Examples of the parametric variations include variations in doping profiles, materials, interconnect width (W ) and thickness (T ), device characteristics like the effective channel length (L e f f ) etc. Physical parameters like W , T of the interconnects and L e f f of the MOS devices vary significantly [START_REF] Boning | Models of process variations in device and interconnect[END_REF] with strong intradie (across die) and inter-die (across wafer) components. These physical variations lead to substantial variations in the electrical parameters viz., conductance, capacitance, inductance, threshold voltages, leakage currents etc., of the CMOS devices and interconnects.

An important aspect of the IC design process is the integrity of the power grid. The exponential increases in transistor density has resulted in huge power distribution networks carrying large transient currents, that result in significant voltage drops (IR and L di dt ) in the voltage levels at the power grid nodes. This combined with continuous reduction in supply voltages (≤ 1.2V for the 90 nm process) make the voltage drops critical as they can severely impact the functionality and the performance of the ICs. Process variations in the power grid conductors and onchip CMOS devices of the functional blocks can have a significant detrimental impact on the voltage levels at the power grid nodes. Given the critical dependence of the sub 100nm ICs on the V DD voltage variations, its important to characterize the effects of process variations on the voltage levels at the power grid nodes.

An established body of literature is available on the analysis and optimization aspects of power grids [1-5, 7, 9-13]. Some of the early work on power grid analysis used DC analysis to obtain the IR drops [START_REF] Dharchoudhury | Design and analysis of power distribution networks in powerPC microprocessors[END_REF][START_REF] Chowdhry | Estimation of maximum currents in MOS IC logic circuits[END_REF][START_REF] Stark | Techniques for calculating currents and voltages in VLSI power supply networks[END_REF]. The variations in current profiles of the functional blocks are often obtained by simulation [START_REF] Dharchoudhury | Design and analysis of power distribution networks in powerPC microprocessors[END_REF][START_REF] Kriplani | Pattern independent minimum current estimation in power and ground buses of CMOS VLSI Circuits[END_REF][START_REF] Krstic | Vector Generation for maximum instantaneous current through supply lines for CMOS circuits[END_REF][START_REF] Chaudhry | Current signature compression for IR-drop analysis[END_REF]. The authors of [START_REF] Nassif | Fast power grid simulation[END_REF] describe a multi-grid approach where coarser meshes are first solved using fast PDE solvers and the results are extrapolated to solve the original mesh. An efficient hierarchical power grid analy-sis technique employing sparsification based on integer linear programming is described in [START_REF] Zhao | Hierarchical analysis of power distribution networks[END_REF] and [START_REF] Qian | Random walks in a supply network[END_REF] presents a linear time algorithm based on random walks that is efficient for incremental localized analysis. The analysis of power grids has only been recently expanded to account for process variations. Assuming normal and lognormal distributions for the device threshold voltages and leakage currents, the authors of [START_REF] Ferzli | Statistical estimation of leakage-induced power grid voltage drop considering within-die process variations[END_REF][START_REF] Ferzli | Statistical Verification of Power Grids Considering Process-Induced Leakage Current Variations[END_REF], present a procedure to compute the mean and bounds on the variance of the voltage drops on a power grid. In [START_REF] Pant | A stochastic approach To power grid analysis[END_REF], the voltage response of a power grid is expressed as a convolution of the grid impulse response and load currents. Their approach to account for variability is to view the load currents of the functional blocks (inputs to the power grid system) as random variables due to the large space of input patterns. Based on this, they compute the mean and variance of the voltage response of the grid.

Our Contributions

In this paper, we propose a new approach [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF] to account for the impact of process variations in the analysis of power grids. Until recently, stochastic analysis of systems generally meant that the system inputs were stochastic but the system itself was deterministic with deterministic parameters. The approach to be described here addresses the important case when the system parameters are also stochastic quantities. For such a case, the physical system itself is an outcome of a stochastic process. Our analysis, thus is completely different from all the previous power grid analyses. In this work, due to the manufacturing variations in the interconnect width W , thickness T and the device channel length L e f f , the electrical parameters in the power grid (R,L,C) are modeled as continuous parameter (spatial) stochastic processes. The key development here is an expansion for the stochastic voltage response as an infinite series of orthogonal polynomials of random variables in an infinite dimensional Hilbert space. The expansion can be optimally truncated to any order depending on the available computational resources and accuracy requirements. With an explicit analytical representation of the voltage response available in terms of the random variables, moments and probability density functions of voltage can be directly computed. This provides an attractive alternative to the computationally expensive Monte Carlo simulations. Our approach has been implemented in a prototype software called OPERA.

In Section 2, we discuss the core contributions of our work. Section 3 contains the problem definition. In Section 4, we discuss the theoretical foundations of our approach. In Section 5, we demonstrate our method through an example followed by the discussion of a special case. We also provide a brief discussion on the implementation issues. The experimental results are shown in Section 6 and conclusions are presented in Section 7.

Problem Statement

We consider an RC model of the power grid with a mesh structure in which the metal interconnects and the vias are modeled as passive RC networks. Power sources provide connections from the external supply to the grid and are modeled as ideal voltage sources between the grid nodes and the ground. The functional blocks (logic gates/latches) that are distributed across the chip act as power drains. They are modeled as known transient current sources between the grid nodes and the ground, in parallel with the non-switching capacitances of the functional blocks. The transient current source profiles are obtained from simulation of the functional blocks at a full supply voltage for a large sequence of random input vectors. The package pin contacts that provide the power supply connection are modeled as resistances in series with the supply sources.

The power grid model can be described in the Laplace domain by the MNA equation

(G + sC) x(s) = U(s) (1)
where (G + sC) is the coefficient matrix and x(s) is the response to be determined. U(s) = (i(s), G 1 V DD ) T is the known excitation; vector i(s) represents the functional block drain currents; G 1 is a diagonal matrix that consists of non-zero elements at those nodes where V DD sources are connected. Almost all the conductance of the power grid comes from the metal interconnects while they contribute only about 5% of the grid capacitance [START_REF] Chaudhry | Design and Analysis of Power Distribution Networks with accurate RLC models[END_REF]. The majority of the capacitance contribution of the power grid comes from the non-switching load capacitances of the gates in the functional blocks. The load capacitance of the gates in turn comes from two major sourcesthe gate (C GS and C GD ) capacitance of the gates they drive and the source/drain diffusion capacitances (C DB and C SB ); the two sources contribute almost equally to the grid capacitance in the present technologies [START_REF] Chaudhry | Design and Analysis of Power Distribution Networks with accurate RLC models[END_REF][START_REF] Rabaey | Digital Integrated Circuits: A Design Perspective[END_REF]. The drain currents i(s) consist of MOS drain currents and the gate oxide and subthreshold leakage currents.

The circuit parameters G and C depend on the grid interconnect and device parameters such as the metal thickness (T ), metal width (W ), channel length (L e f f ) etc.. In addition, the MOS drain currents i(s) in the excitation U(s) vary significantly with changes in L e f f . In the presence of random process variations, these physical characteristics of the interconnects and devices, and hence their electrical characteristics are modeled as spatial stochastic processes (they vary randomly and spatially across a die for intra-die variations and across a wafer for inter-die variations).

In the present work, we assume that the conductance G of the grid varies with parameters W and T . We ignore variations in C due to W and T . We assume that the gate capacitance contributes about 40% of the total grid capacitance and that this gate capacitance varies with L e f f . We have from SPICE models C gate ∝ W e f f L e f f C ox where W e f f ,C ox represent the effective width and gate oxide capacitance of a MOS transistor. The drain currents i(s) that consist of MOS drain currents and leakage currents are known to vary significantly with L e f f . We consider only the inter-die variations in this work and hence the variations in G and C are fixed for any one single die. Thus, they are modeled as random variables and not as spatial random processes.

Consider the manufacturing process. Every trial of the manufacturing process (trial denotes the fabrication of a die or a wafer) results in a different value for W , T and L e f f of the interconnects and devices. These parameters are thus functions that map one point in the manufacturing sample space to some real value. Thus, W , T and L e f f are random variables over the manufacturing sample space. In general, let Ω denote the manufacturing sample space. For ω ∈ Ω, ξ i : ω → R denotes a random variable. Let ξ ξ ξ(ω) = (ξ 1 (ω),...,ξ n (ω)) denote a vector of n such random vari-ables. Let Θ : Ω → R denote the vector space of all random variables ξ i . In the presence of process variations the MNA equations for the interconnect can be expressed as:

( G(ξ ξ ξ(ω)) + s C(ξ ξ ξ(ω)) ) x(s, ξ ξ ξ(ω)) = U(s, ξ ξ ξ(ω)) (2)
Equation ( 2) is known as a stochastic differential equation as the operator (G(ξ ξ ξ(ω)) + sC(ξ ξ ξ(ω))) is a stochastic process dependent on a random variable ξ ξ ξ and the deterministic Laplace parameter s. Further, the excitation vector U(s, ξ ξ ξ(ω)) has both deterministic and random components. For each manufacturing outcome ω and for each corresponding value of the parameter ξ(ω), x(s, ξ ξ ξ(ω)) denotes the fixed response of the system for that particular manufacturing outcome. In the next section, we discuss the theoretical foundations of our proposed approach to compute the response x(s, ξ ξ ξ(ω)).

Proposed Approach

The approach presented here is based on representing the stochastic voltage response x(s, ξ ξ ξ(ω)) of the power grid as an infinite series of orthogonal polynomials in an infinite dimensional Hilbert space of random variables. We assume that x(s, ξ ξ ξ(ω)) is a second order process, i.e. all the random variables have finite variances. Without loss of generality, we also assume that the variables ξ ξ ξ(ω) are random variables with zero mean and unit variance. Given a random variable W with mean W µ and standard deviation W σ , W can always be expressed in a normalized form as W = W µ + W σ ξ W , where ξ W is a random variable with zero mean and unit variance.

For simplicity, henceforth we will write ξ ξ ξ(ω) as ξ ξ ξ. To explain our approach, we state a few well known facts from the theory of orthonormal expansions [START_REF] Wong | Stochastic Processes in Information and Dynamical Systems[END_REF].

• Let V be an inner product space, with the inner product denoted by •, • . For non-zero x, y ∈ V , x and y are orthogonal if x, y = 0. They are orthonormal if x = y = 1. • An complete inner product space H (i.e. one in which every Cauchy sequence converges) is called a Hilbert space. • An infinite family of orthonormal vectors, {φ} ∞ k=1 , in a Hilbert space is called an orthonormal basis if it is a maximal set of mutually orthonormal vectors.

• If {φ} ∞ k=1 is a orthonormal basis of a Hilbert space H, then the infinite series ∑ ∞ k=1 x, φ k φ k converges in norm to x. We now return to the problem of representing the stochastic voltage response x(s, ξ ξ ξ). From Section 3, we know that the space Θ : {Ω → R} denotes the infinite dimensional vector space of mappings, each mapping representing a random variable ξ ξ ξ(ω). Let P be the probability measure on the sample space Ω i.e., the random variables ξ ξ ξ have the probability density function P. This space of mappings forms a Hilbert space H Θ , where the inner product of any two mappings is the expected value of their product under P. That is,

< ξ ξ ξ m , ξ ξ ξ n >= E(ξ ξ ξ m , ξ ξ ξ n ) = Z Ω ξ ξ ξ m ξ ξ ξ n dP (3)
The stochastic voltage response x(s, ξ ξ ξ) is an element of Θ. Therefore, if we can find a orthonormal basis {γ 1 (ξ ξ ξ), γ 2 (ξ ξ ξ),...} for H Θ , then x(s, ξ ξ ξ) can be represented by

x(s, ξ ξ ξ) = ∞ ∑ n=0 a n (s) γ n (ξ ξ ξ) (4)
Equation ( 4) is a general representation of the stochastic voltage response process x(s, ξ ξ ξ) in terms of the orthonormal basis functions of random variables ξ ξ ξ. The primary task now is to identify the orthonormal basis {γ 1 (ξ ξ ξ), γ 2 (ξ ξ ξ),...} in ξ ξ ξ. One such orthonormal basis is the set of Hermite polynomials of all orders. This basis is valid for any second order process. However, other types of polynomials can also serve as orthonormal bases. The orthogonality is defined w.r.t the norm in Equation (3) which is dependent on the probability density function P of ξ ξ ξ. So, for different probability distributions of the random variables, different orthonormal basis sets need to be identified. The well established Askey scheme of polynomials [START_REF] Xiu | Modeling uncertainty in flow simulations via generalized polynomial chaos[END_REF] helps us identify the orthogonal polynomials for different probability density functions. For example, if the underlying random variables are Gaussian or lognormal, then the best choice (in terms of speed of convergence) would be Hermite polynomials. Similarly, for Gamma, Beta and Uniform random variables, the best choices would be Laguerre, Jacobi and Legendre polynomials respectively.

For the sake of demonstrating our approach, we consider a power grid subjected to inter-die Gaussian variations in interconnect width (W ), thickness (T ) and device channel length (L e f f ). The methodology described here can be easily extended to consider other probability distributions and any number of variables.

Gaussian Random Variables

Let ξ ξ ξ = {ξ 1 , ξ 2 , ξ 3 ,...,ξ n } ∈ Θ be zero mean orthonormal Gaussian random variables. Hermite polynomials of all orders in ξ ξ ξ form an orthonormal basis for H Θ . They are defined by

H p ({α 1 , α 2 , ••• , α p }) = (-1) p e 1 2 ξ ξ ξ t ξ ξ ξ ∂ p ∂α 1 ∂α 2 •••∂α p e -1 2 ξ ξ ξ t ξ ξ ξ , (5)
where {α i } are any set of p variables chosen from the set {ξ 1 , ξ 2 , ••• , ξ n } with repetitions. Since p variables can be chosen from a set of n variables in M = (p + n -1)!/p!(n -1)! ways, the number of Hermite polynomials of degree p is given by M. As an example, Hermite polynomials in {ξ 1 , ξ 2 } of order 0, 1, 2 can be found as follows.

order 0: H 0 ({}) = 1, order 1:

H 1 (ξ 1 ) = ξ 1 , H 1 (ξ 2 ) = ξ 2 , order 2: H 2 (ξ 1 , ξ 1 ) = ξ 2 1 -1, H 2 (ξ 1 , ξ 2 ) = ξ 1 ξ 2 , H 2 (ξ 2 , ξ 2 ) = ξ 2 2 -1 (6)
Let {ξ 1 , ξ 2 , ξ 3 ,...} denote an infinite set of zero mean orthonormal Gaussian random variables. Then the stochastic response can be represented as a convergent infinite series expansion of Hermite polynomials as [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF] x(s,

ξ ξ ξ) = b 0 (s) H 0 + ∞ ∑ i 1 =1 b i 1 (s) H 1 (ξ i 1 ) + ∞ ∑ i 1 =1 i 1 ∑ i 2 =1 b i 1 i 2 (s)H 2 (ξ i 1 , ξ i 2 ) ( 7
)
+ ∞ ∑ i 1 =1 i 1 ∑ i 2 =1 i 2 ∑ i 3 =1 b i 1 i 2 i 3 (s)H 3 (ξ i 1 , ξ i 2 , ξ i 3 )
. . .

Note:

A closer look at Equations ( 6), ( 7), [START_REF] Nassif | Fast power grid simulation[END_REF] shows that Equation ( 4) is simply a rearrangement of the terms in Equation ( 7). The remaining task then is to determine the coefficients {a i } of the expansion given in Equation ( 4). The expansion shown in Equation ( 7) called the Homogeneous Chaos was developed by Weiner [START_REF] Wiener | The Homogeneous Chaos[END_REF]. The celebrated result of Cameron and Martin [START_REF] Cameron | The orthogonal development of non-linear functionals in series of Fourier-Hermite functionals[END_REF] extended the expansions to general function spaces, and is now referred to as Polynomial Chaos. Ghanem and Spanos [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF] developed applications of these results to the study of systems with stochastic parameters.

Evaluation of the Coefficients

Though the number of random variables ξ ξ ξ are finite, the order p and hence the number of terms N + 1 in the stochastic response expansion are infinite. For practical computation purposes, the infinite series needs to be projected onto a finite space, i.e. truncated. The general approach is to limit the response expansion to a finite order p, which in turn determines the accuracy. Often, a second order (p = 2) or third order (p = 3) expansion is sufficient. If there are n random variables, the response expansion obtained by limiting the order to p would be

x ( s, ξ ξ ξ) = N ∑ i=0 a i (s) γ i (ξ ξ ξ) (8)
where

N = ∑ p k=0 n -1 + k k .
The error due to truncation is given by

∆ p (s, ξ ξ ξ) = (G(ξ ξ ξ) + s C(ξ ξ ξ)) x(s, ξ ξ ξ) -U(s, ξ ξ ξ) (9)
Once we have the truncated expansion from Equation ( 8), we need to evaluate the best deterministic coefficients {a i } that result in the best minimization of the truncation error. We follow the principle of orthogonality which states that the best approximation of the response x(s, ξ ξ ξ) is one in which the truncation error ∆ p (s, ξ ξ ξ) is orthogonal to the approximation. The application of the principle of orthogonality to obtain a finite projection is known as the Galerkin method. The orthogonality is defined w.r.t to some norm, which in this case is given by < ∆ p (s, ξ ξ ξ), γ i (ξ ξ ξ) >. Thus, to obtain the deterministic coefficients we need to solve the system of equations given by < ∆ p (s, ξ ξ ξ), γ i (ξ ξ ξ) >= 0, i = 0, 1, 2, •••N (10)

Illustration

As an example to illustrate the methods, we consider a power grid subjected to random process variations where ξ W , ξ T , ξ L denote the normalized variations in width, thickness and the effective device channel length respectively. Then

ξ W = ∆W W , ξ T = ∆T T , ξ L = ∆L e f f L e f f ( 11
)
Without loss of generality, we assume that they are uncorrelated Gaussian random variables. Certainly if they were not, given their covariance matrix, they can always be transformed into a set of uncorrelated random variables by an orthogonal transformation technique like principal component analysis [START_REF] Papoulis | Probability, Random Variables and Stochastic Processes[END_REF]. As discussed earlier, the G matrix depends on ξ W and ξ T and the C matrix depends on ξ L . We use a linear model to capture the dependence of G and C on these random variables in accordance with the models in contemporary literature [START_REF] Mehrotra | A methodology for modeling the effects of systematic within-die interconnect and device variation on circuit performance[END_REF]. But there are no limitations on the specific model to be chosen. The MNA equation for the grid system is given by:

(G(ξ ξ ξ) + sC(ξ ξ ξ)) x(s, ξ ξ ξ) = U(s, ξ ξ ξ) (12)
where ξ ξ ξ = (ξ W , ξ T , ξ L ) is the random variable vector and the perturbed matrices G and C and U(s) are given by

G(ξ ξ ξ) = G a (s) + G b (s) ξ W + G c (s) ξ T C(ξ ξ ξ) = C a (s) +C c (s) ξ L (13) U(s, ξ ξ ξ) = U a (s) +U b (s) ξ W +U d (s) ξ T +U c (s) ξ L G a ,C a ,U a are the mean matrices and G b , G c , C c , U b , U c , U d repre- sent the perturbation matrices of G, C, U(s) w.r.t W , T , L e f f .
Consider Equation ( 13) and the physical definition for conductance per unit length G = W T ρ . We can observe that for the linear model in Equation [START_REF] Ferzli | Statistical Verification of Power Grids Considering Process-Induced Leakage Current Variations[END_REF], G b and G c are same as G a scaled by some constants. So, we have G b = d G a and G c = eG a where e and d are some constants. Since the scaled sum of two independent Gaussian variables ξ W and ξ T (d ξ W + e ξ T ) is a Gaussian variable with calcuable mean and variance, matrix G and hence U(s, ξ ξ ξ) can be re-represented using a single normalized variable ξ G as

G(ξ ξ ξ) = G a (s) + G g (s) ξ G , U(ξ ξ ξ) = U a (s) +U g (s) ξ G +U c (s) ξ L (14) C(ξ ξ ξ) = C a (s) +C c (s) ξ L ξ ξ ξ is now given by ξ ξ ξ = (ξ G , ξ L ) T .
Using Hermite polynomials as the basis, we can expand the response x(s, ξ ξ ξ) using a second order expansion (p = 2) from Equation ( 7) as

x(s, ξ ξ ξ) = a 0 (s) + a 1 (s) ξ G + a 2 (s) ξ L + a 3 (s)(ξ G 2 -1) + a 4 (s)(ξ G ξ L ) + a 5 (s)(ξ L 2 -1). ( 15
)
where a i (s) is a vector. The coefficient vectors {a i } need to be determined.

Following the orthogonal truncation method from Section 4.2, we can define the error as in [START_REF] Krstic | Vector Generation for maximum instantaneous current through supply lines for CMOS circuits[END_REF].

∆ p (s, ξ ξ ξ) = (G(ξ ξ ξ) + sC(ξ ξ ξ)) x(s, ξ ξ ξ) -U(s, ξ ξ ξ) (16)
Coefficient vectors {a i (s)} are obtained by solving (see Equation [START_REF] Chaudhry | Current signature compression for IR-drop analysis[END_REF])

∆ p (s, ξ ξ ξ), γ j (ξ ξ ξ) = 0 for j = 0, 1, 2,...,N (17)
The inner product Σ p (s, ξ ξ ξ), γ j is defined as

∆ p (s, ξ ξ ξ), γ j (ξ ξ ξ) = Z +∞ -∞ Z +∞ -∞ ∆ p (s, ξ ξ ξ)γ j (ξ ξ ξ)W (ξ ξ ξ) dξ G dξ L = 0, ( 18
)
where W (ξ ξ ξ) is the bivariate Gaussian probability density function.

Thus for each j = 0, 1, 2, ••• , N, Equation ( 17) results in a linear system of (N + 1) equations to solve for the deterministic coefficients vectors {a i } represented by the vector a(s).

( G + s C) a(s) = Ũ(s), (19)
where

G =         G a G g 0 0 0 0 G g G a 0 2G g 0 0 0 0 G a 0 G g 0 0 2G g 0 2G a 0 0 0 0 G g 0 G a 0 0 0 0 0 0 2G a         (20)
C =         C a 0 C c 0 0 0 0 C a 0 2C b C c 0 C c 0 C a 0 0 2C c 0 0 0 2C a 0 0 0 C c 0 0 C a 0 0 0 2C c 0 0 2C a         (21)
Matrix Ũ(s) is given by

Ũ(s) = (U a (s), U g (s), U c (s), 0, 0, 0 ) T (22)
Now, we can solve Equation [START_REF] Xiu | Modeling uncertainty in flow simulations via generalized polynomial chaos[END_REF] numerically to obtain the coefficient vector a(s). Once the vector a(s) is obtained, we have an explicit expression for the circuit response x(s, ξ ξ ξ) in terms of ξ ξ ξ given by Equation [START_REF] Rabaey | Digital Integrated Circuits: A Design Perspective[END_REF]. With this explicit expression, we can obtain the moments of the stochastic voltage response at any node of the power grid as follows:

Mean(x(t, ξ ξ ξ)) = a 0 (t) Var(x(t, ξ ξ ξ)) = (a 1 (t)) 2 Var(ξ G ) + (a 2 (t)) 2 Var(ξ L ) + (a 3 (t)) 2 Var(ξ G 2 -1) + (a 4 (t)) 2 Var(ξ G ) Var(ξ L ) + (a 5 (t)) 2 Var(ξ L 2 -1) (23) Var(x(t, ξ ξ ξ)) = (a 1 (t)) 2 + (a 2 (t)) 2 + (a 3 (t)) 2 2 + (a 4 (t)) 2 + (a 5 (t)) 2 2
To obtain higher order moments for x(t, ξ ξ ξ), we can use the equality E(x n (t, ξ ξ ξ)) = < x n-1 (t, ξ ξ ξ), x(t, ξ ξ ξ) >, provided that x(t, ξ ξ ξ) has an accurate representation using an expansion of sufficiently large order p. Once the higher order moments are obtained, expansions like Gram-Charlier series or Edgeworth series could be used to obtain the probability density function of x(t, ξ ξ ξ) directly.

Special Case

As a special case, let's consider only the variations in the drain currents U(s), i.e. the R.H.S of the MNA analysis Equation (1) say due to threshold variations. U(s) consists of MOS drain currents and leakage currents; the latter are known to vary exponentially with the threshold voltages. Our problem then becomes similar to the one addressed by Ferzli and Najm [START_REF] Ferzli | Statistical estimation of leakage-induced power grid voltage drop considering within-die process variations[END_REF][START_REF] Ferzli | Statistical Verification of Power Grids Considering Process-Induced Leakage Current Variations[END_REF]. We then have,

(G + sC) x(s, ξ ξ ξ) = U(s, ξ ξ ξ) (24)
where ξ ξ ξ represents the threshold variations. If a normal distribution for V th is assumed, then the distribution for the leakage currents becomes lognormal. To consider intra-die variations in V th , let's divide the chip in to a finite number of regions say 2 for this example. Let's assume that ξ ξ ξ = (ξ 1 , ξ 2 ) are the normalized uncorrelated Gaussian random variables that represent V th variations for the two chip regions. Then, the currents U(s, ξ ξ ξ) can always be expressed using an orthogonal Hermite polynomial basis [START_REF] Debusschere | Numerical challenges in the use of Polynomial Chaos representations for stochastic processes[END_REF] in ξ 1 , ξ 2 to any required order (p) of accuracy, say p = 2 for this case. Since, the R.H.S of Equation ( 24) is stochastic, the response of the grid becomes stochastic as well and can be represented using the Hermite basis. We thus have,

x(s, ξ ξ ξ) = x 0 (s) + x 1 (s) ξ 1 + x 2 (s) ξ 2 + x 3 (s)(ξ 1 2 -1)

+ x 4 (s)(ξ 1 ξ 2 ) + x 5 (s)(ξ 2 2 -1) (25) U(s, ξ ξ ξ) = U 0 (s) +U 1 (s) ξ 1 +U 2 (s) ξ 2 +U 3 (s)(ξ 1 2 -1) + U 4 (s)(ξ 1 ξ 2 ) +U 5 (s)(ξ 2 2 -1) (26)
Following the error minimization procedure illustrated in Section 5, our analysis simply translates to solving independent equations of the form for n = 0, 1,...,5

(G + sC) x n (s) = U n (s) (27)
All we need is a single LU factorization of the original matrix (G + sC) and then repeated solves for different values of the R.H.S. Using the formulae from Equation ( 23), we can directly compute the mean, variance and other higher order moments unlike [START_REF] Ferzli | Statistical estimation of leakage-induced power grid voltage drop considering within-die process variations[END_REF][START_REF] Ferzli | Statistical Verification of Power Grids Considering Process-Induced Leakage Current Variations[END_REF] which can calculate only some bounds for the variance.

Implementation Issues

The procedure for obtaining the voltage response has been implemented in a prototype software called OPERA. One of the primary issues involved is the computational complexity of our approach. Solving the system of Equations ( 19) is computationally the most intensive step in our approach. The complexity of that step is dependent on the length of the vector a(s) which varies as O(r p ), where p is the order of the expansion and r is the number of random variables. We found an order 2/order 3 expansion to be sufficiently accurate for variational power grid analysis considering realistic bounds for maximum variability in the grid interconnect and device parameters. Also, the matrices G and C are very sparse and they have been observed to become increasingly sparser with an increase in the order of the expansion (p) or the number of random variables (r).

Further, computational complexity of OPERA can be significantly reduced by efficient techniques like model order reduction (MOR) [START_REF] Feldmann | Sparse and Efficient Reduced Order Modeling of Linear Subcircuits with Large Number of Terminals[END_REF], multi grid analysis [START_REF] Nassif | Fast power grid simulation[END_REF] and iterative block solvers with appropriate pre-conditioners [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF]. MOR techniques can be used as the power grid node voltages in the top layers and their moments w.r.t ξ ξ ξ are typically of no interest to the designer. Model order related stability issues have been addressed in a number of literatures and any existing stability technique can be incorporated in our method.

Experimental Results

OPERA has been verified for many industrial power grids considering order 2/order 3 expansions for representing the stochastic voltage response. Results for a few grids are presented here. As discussed earlier, we assume a Gaussian distribution for the variables W , T and L e f f for all the grids. We consider a linear or a first order model for the variations in G, C and the drain currents i(s). Leakage currents are known to vary exponentially with L e f f which may suggest a higher order expansion for U(s) in L e f f . But for the purposes of this paper, we consider a linear expansion itself as they constitute only 5% of the total currents in the current CMOS technologies [START_REF]The International Technology Roadmap for Semiconductors[END_REF].

Table 1 shows the results for the transient analysis of 7 industrial grids for maximum 3σ variations of 20% in ξ W , 15% in ξ T (hence 25% in ξ G ) and 20% in ξ L . Note that we are considering the inter-die variations and that only 40% of the capacitance varies with L e f f as discussed earlier. A fixed time step was used in carrying out the transient analysis and an order 2 expansion was used for representing the stochastic response. Comparison of results from OPERA with Monte Carlo simulations (1000 samples for each case) has been done for all the grids and the results are reported in Table 1. The average and the maximum errors in obtaining the mean (µ) and variance (σ) of the voltage response from OPERA compared to Monte Carlo simulations are shown in the table. The errors reported in Table 1 are for data obtained from simulation across all nodes and all time points of the transient simulation of the grid. We can see from the table that OPERA demonstrates good accuracy in determining the mean and variance of the voltage response. The accuracy in obtaining the variance by OPERA can be increased further by increasing the order of the expansion.

The drain current profiles used for the transient analysis of the power grid were such that the peak drop in the voltage at any grid node was less than 10 % of the V DD . Under such conditions, it was observed that for all grids the mean voltage drops at the grid nodes (µ) with variations was more or less the same as the nominal voltage drops (µ 0 ) with out variations. And the difference between the two drops (µµ 0 ) was negligible when expressed as a fraction of % of V DD . However on average for each grid, the ± 3σ variation in the voltage drops at the grid nodes, was about ± 35 % of µ 0 , where µ 0 are the nominal voltage drops with out variations. This strongly supports the necessity of considering the effects of process variations on power grids.

For the power grid with 19,181 nodes, we plotted the distribution of the voltage response from OPERA and Monte Carlo simulations w.r.t to the variations in ξ W , ξ T and ξ L at arbitrarily selected nodes in the power grid. Figures 1,2 show the voltage distribution at a select node from Monte Carlo simulations and OPERA.

Conclusions

We presented a general framework to accurately compute the stochastic response of power grids in the presence of process variations. Using orthogonal polynomial expansions in a Hilbert space, we provide an explicit analytical representation of the stochastic response. The coefficients of the analytical expansion are determined precisely by utilising the orthogonality property of the polynomials of the expansion. Further, the expansion facilitates the direct and precise computation of the moments of the power grid response. We also show that our analysis becomes very simplified if we consider just the impact of leakage current variations on the power grid. We implemented the algorithm in a prototype software called OPERA and verified it extensively for many industrial grids. OPERA demonstrates very good accuracy when compared with the classical Monte Carlo simulations along with providing significant speed-ups up to two orders of magnitude.

Figure 1 .

 1 Figure 1. Voltage distribution from OPERA and MC

Figure 2 .

 2 Figure 2. Voltage distribution from OPERA and MC

Table 1 . Results for grids from OPERA and Monte Carlo simulations for order 2 expansion

 1 SizeAve. % Error Max. % Error Ave. % Error Max. % Error ± 3σ variation

	CPU time	CPU time	Speedup

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

Acknowledgements

This work was carried out under the auspices of the National Science Foundation's Grant #CCR-0205227 and Grant #EEC-9523338. It was conducted at the National Science Foundation's State/Industry/University Cooperative Research Centers' (NSF-S/IUCRC) Center for Low Power Electronics (CLPE) which is supported by the NSF, the State of Arizona, and an industrial consortium. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.