Enabling High-Dimensional Hierarchical Uncertainty Quantification by ANOVA and Tensor-Train

14 Pages·2014·1.47 MB·English

IEEETRANSACTIONSONCOMPUTER-AIDEDDESIGNOFINTEGRATEDCIRCUITSANDSYSTEMS, VOL. XX,NO.XX,XX2015 1
Enabling High-Dimensional Hierarchical
Uncertainty Quantification by ANOVA and
Tensor-Train Decomposition
Zheng Zhang, Xiu Yang, Ivan V. Oseledets, George Em Karniadakis, and Luca Daniel
Accepted by IEEE Trans. Computer-Aided Design of Integrated Circuits and Systems
Abstract—Hierarchical uncertainty quantification can reduce fast convergence rate, such techniques have been successfully
the computational cost of stochastic circuit simulation by em- applied in the stochastic analysis of integrated circuits [14]–
ploying spectral methods at different levels. This paper presents
[20], VLSI interconnects [21]–[25], electromagnetic [26] and
an efficient framework to simulate hierarchically some chal-
MEMS devices [1], [27], achieving significant speedup over
lenging stochastic circuits/systems that include high-dimensional
subsystems. Due to the high parameter dimensionality, it is Monte Carlo when the parameter dimensionality is small or
challenging to both extract surrogate models at the low level medium.
of the design hierarchy and to handle them in the high-level Sincemanyelectronicsystemsaredesignedinahierarchical
simulation. In this paper, we develop an efficient ANOVA-
way,itispossibletoexploitsuchstructureandsimulateacom-
based stochastic circuit/MEMS simulator to extract efficiently
plex circuit by hierarchical uncertainty quantification [28]1.
thesurrogatemodelsatthelowlevel.Inordertoavoidthecurse
of dimensionality, we employ tensor-train decomposition at the Specifically, one can first utilize stochastic spectral methods
highleveltoconstructthebasisfunctionsandGaussquadrature to extract surrogate models for each block. Then, circuit
points. As a demonstration, we verify our algorithm on a equations describing the interconnection of blocks may be
stochasticoscillatorwithfourMEMScapacitorsand184random
solvedwithstochasticspectralmethodsbytreatingeachblock
parameters. This challenging example is simulated efficiently by
our simulator at the cost of only 10 minutes in MATLAB on a as a single random parameter. Typical application examples
regular personal computer. include (but are not limited to) analog/mixed-signal systems
(e.g., phase-lock loops) and MEMS/IC co-design. In our pre-
Index Terms—Uncertainty quantification, hierarchical uncer-
liminary conference paper [1], this method was employed to
tainty quantification, generalized polynomial chaos, stochastic
modeling and simulation, circuit simulation, MEMS simulation, simulatealow-dimensionalstochasticoscillatorwith9random
highdimensionality,analysisofvariance(ANOVA),tensortrain. parameters, achieving 250 speedup over the hierarchical
×
Monte-Carlo method proposed in [32].
Paper Contributions. This paper extends the recently
I. INTRODUCTION developed hierarchical uncertainty quantification method [28]
PROCESS variations have become a major concern in to the challenging cases that include subsystems with high
dimensionality (i.e., with a large number of parameters). Due
submicron and nano-scale chip design [2]–[6]. In or-
to such high dimensionality, it is too expensive to extract a
der to improve chip performances, it is highly desirable to
surrogatemodelforeachsubsystembyanystandardstochastic
develop efficient stochastic simulators to quantify the un-
spectral method. It is also non-trivial to perform high-level
certainties of integrated circuits and microelectromechanical
simulation with a stochastic spectral method, due to the high-
systems (MEMS). Recently, stochastic spectral methods [7]–
dimensional integration involved when computing the basis
[12] have emerged as a promising alternative to Monte Carlo
functions and Gauss quadrature rules for each subsystem. In
techniques [13]. The key idea is to represent the stochastic
order to reduce the computational cost, this work develops
solution as a linear combination of some basis functions
some fast numerical algorithms to accelerate the simulations
(e.g.,generalizedpolynomialchaos[8]),andthencomputethe
at both levels:
solutionbystochasticGalerkin[7],stochasticcollocation[9]–
[12] or stochastic testing [14]–[16] methods. Due to the At the low level, we develop a sparse stochastic testing
•
simulator based on adaptive anchored ANOVA [33]–[37]
Somepreliminaryresultsofthisworkhavebeenreportedin[1].Thiswork
to efficiently simulate each subsystem. This approach
was funded by the MIT-SkolTech program. I. Oseledets was also supported
bytheRussianScienceFoundationunderGrant14-11-00659. exploits the sparsity on-the-fly, and it turns out to be
Z. Zhang and L. Daniel are with the Research Laboratory of Electronics, suitable for many circuit and MEMS problems. This
Massachusetts Institute of Technology (MIT), Cambridge, MA 02139, USA
algorithm was reported in our preliminary conference
(e-mail:z [email protected],[email protected]).
X.YangwaswiththeDivisionofAppliedMathematics,BrownUniversity, paper [1] and was used for the global sensitivity analysis
Providence, RI 02912. Now he is with the Pacific Northwest National of analog integrated circuits.
Laboratory,Richland,WA99352,USA(e-mail:[email protected]).
G.KarniadakisiswiththeDivisionofAppliedMathematics,BrownUniver-
sity,Providence,RI02912,USA(e-mail:george [email protected]). 1Design hierarchy can be found in many engineering fields. In the recent
IvanV.OseledetsiswiththeSkolkovoInstituteofScienceandTechnology, work[29]ahierarchicalstochasticanalysisandoptimizationframeworkbased
Skolkovo143025,Russia(e-mail:[email protected]). onmulti-fidelitymodels[30],[31]wasproposedforaircraftdesign.
r
In the high-level stochastic simulation, we accelerate the h
•
three-termrecurrencerelation[38]bytensor-traindecom-
position[39]–[41].Ouralgorithmhasalinearcomplexity Higher-Level Eq.
with respect to the parameter dimensionality, generating
y y
a set of basis functions and Gauss quadrature points 1 y q
surrogate i
withhighaccuracy(closetothemachineprecision).This surrogate
model ... surrogate ...
algorithm was not reported in [1]. r model r rmodel
x x x
1 i q
II. BACKGROUNDREVIEW lower-level lower-level lower-level
parameters parameters parameters
This section first reviews the recently developed stochastic
testing circuit/MEMS simulator [14]–[16] and hierarchical
Fig.1. Demonstrationofhierarchicaluncertaintyquantification.
uncertaintyquantification[28].Thenweintroducesomeback-
ground about tensor and tensor decomposition.
In order to compute ~x(t,ξ~), stochastic testing [14]–[16]
A. Stochastic Testing Circuit/MEMS Simulator substitutes x˜(t,ξ~) into (1) and forces the residual to zero at
K testing samples of ξ~. This gives a deterministic differential
Given a circuit netlist (or a MEMS 3D schematic file),
algebraic equation of size nK
device models and process variation descriptions, one can set
up a stochastic differential algebraic equation: dq(ˆx(t))
+f(ˆx(t),u(t))=0, (6)
dt
d~q ~x(t,ξ~),ξ~
(cid:16) (cid:17) +f~ ~x(t,ξ~),ξ~,u(t) =0 (1) where the state vector ˆx(t) contains all coefficients in (3).
dt (cid:16) (cid:17) Stochastic testingthen solves Eq.(6)withalinearcomplexity
where ~u(t) is the input signal, ξ~=[ξ , ,ξ ] Ω Rd are ofK andwithadaptivetimestepping,andithasshownhigher
1 d
d mutually independent random variab·l·e·s desc∈ribin⊆g process efficiency than standard stochastic Galerkin and stochastic
variations. The joint probability density function of ξ~ is collocation methods in circuit simulation [14], [15].
1) ConstructingBasisFunctions: ThebasisfunctionH (ξ~)
α~
d
is constructed as follows (see Section II of [16] for details):
ρ(ξ~)= ρ (ξ ), (2)
k k
First,forξ oneconstructsasetofdegree-α orthonormal
kY=1 • univariateipolynomials ϕi (ξ ) p acciording to its
wξhereΩρk.(ξInk)ciisrctuhiet manaarlgyisniasl, p~xrobRanbildietynodteenssnitoydfaulnvcotilotangeosf marginal probability den(cid:8)sityαiρi(iξi(cid:9))α.i=0
aknd∈brankch currents; ~q Rn and∈f~ Rn represent charge/flux • Next, based on the obtained univariate polynomials of
∈ ∈ each random parameter one constructs the multivariate
andcurrent/voltage,respectively.InMEMSanalysis,Eq.(1)is basis function: H (ξ~)= d ϕi (ξ ).
theequivalentformofacommonlyused2nd-orderdifferential α~ i=1 αi i
The obtained basis functioQns are orthonormal polynomials
equation [1], [42]; ~x includes displacements, rotations and
in the multi-dimensional parameter space Ω with the density
their first-order derivatives with respect to time t.
When~x(ξ~,t)hasaboundedvarianceandsmoothlydepends measure ρ(ξ~). As a result, some statistical information can be
on ξ~, we can approximate it by a truncated generalized easily obtained. For example, the mean value and variance of
~x(t,ξ~) are xˆ (t) and (xˆ (t))2, respectively.
polynomial chaos expansion [8] 0 α~
2) Testing Point Seα~lPe6=c0tion: The selection of testing points
~x(t,ξ~) x˜(t,ξ~)= xˆ (t)H (ξ~) (3)
≈ α~ α~ influencethenumericalaccuracyofthesimulator.Instochastic
α~X∈P testing,thetestingpoints ξ~j K areselectedbythefollowing
where xˆα~(t) Rn denotes a coefficient indexed by vector two steps (see Section III{-C}ojf=[114] for details):
α~ = [α , ,∈α ] Nd, and the basis function H (ξ~) is a
1 d α~ First, compute a set of multi-dimensional quadrature
··· ∈
multivariate polynomial with the highest order of ξ being α . •
i i points. Such quadrature points should give accurate re-
In practical implementations, it is popular to set the highest
sultsforevaluatingthenumericalintegrationofanymul-
total degree of the polynomials as p, then = α~ αk tivariate polynomial of ξ~ over Ω [with density measure
N, 0 α1 + +αd p and the total Pnumbe{r o|f basi∈s ρ(ξ~)] when the polynomial degree is 2p.
≤ ··· ≤ }
functions is ≤
Next, among the obtained quadrature points, we select
•
p+d (p+d)! the K samples with the largest quadrature weights under
K = = . (4)
(cid:18) p (cid:19) p!d! the constraint that V RK×K is well-conditioned. The
(j,k) element of V is∈H (ξ~j).
For any integer j in [1,K], there is a one-to-one correspon- k
dence between j and α~. As a result, we can denote a basis
function as H (ξ~) and rewrite (3) as B. Hierarchical Uncertainty Quantification
j
Consider Fig. 1, where an electronic system has q sub-
K
~x(t,ξ~) x˜(t,ξ~)= xˆj(t)Hj(ξ~). (5) systems. The output yi of a subsystem is influenced by
≈ Xj=1 some process variations ξ~i Rdi, and the output ~h of the
∈
whole system depends on all random parameters ξ~’s. For
i
simplicity, in this paper we assume that y only depends
i
on ξ~ and does not change with time or frequency. Directly
i
simulatingthewholesystemcanbeexpensiveduetothelarge
problem size and high parameter dimensionality. If y ’s are
i
mutuallyindependentandsmoothlydependentonξ~’s,wecan
i
accelerate the simulation in a hierarchical way [28]:
First, perform low-level uncertainty quantification. We
•
use stochastic testing to simulate each block, obtaining Fig. 2. Demonstration of a vector (left), a matrix (center) and a 3-mode
a generalized polynomial expansion for each y . In this tensor(right).
i
stepwecanalsoemployotherstochasticspectralmethods
such as stochastic Galerkin or stochastic collocation.
samplesarerequiredtoobtainy .Second,whenhighaccuracy
Next, perform high-Level uncertainty quantification. By i
• is required, it is expensive to implement (7) due to the non-
treating y ’s as the inputs of the high-level equation,
i trivial integrals when computing κ and γ . Since the density
we use stochastic testing again to efficiently compute j j
function of ζ is unknown, the integrals must be evaluated in
~h. Since y has been assumed independent of time and i
i the domain of ξ~, with a cost growing exponentially with d
frequency, we can treat it as a random parameter. i i
when a deterministic quadrature rule is used.
In order to apply stochastic spectral methods at the high
level, we need to compute a set of specialized orthonormal
C. Tensor and Tensor Decomposition
polynomials and Gauss quadrature points/weights for each
inputrandomparameter.Forthesakeofnumericalstability,we Definition 1 (Tensor). A tensor A RN1×N2×···×Nd is a
∈
define a zero-mean unit-variance random variable ζ for each multi-mode (or multi-way) data array. The mode (or way) is
i
subsystem, by shifting and scaling y . The intermediate vari- d, the number of dimensions. The size of the k-th dimension
i
ables ζ~ = [ζ1,··· ,ζq] are used as the random parameters in is Nk. An element of the tensor is A(i1,··· ,id), where the
the high-level equation. Dropping the subscript for simplicity, positiveintegerik istheindexforthek-thdimensionand1
≤
wedenoteageneralintermediate-levelrandomparameterbyζ ik Nk.ThetotalnumberofelementsofAisN1 Nd.
≤ ×···×
and its probability density function by ρ(ζ) (which is actually
As a demonstration, we have shown a vector (1-mode
unknown),thenwecanconstructp+1orthogonalpolynomials tensor) in R3 1, a matrix (2-mode tensor) in R3 3 and a
{πj(ζ)}pj=0 via a three-term recurrence relation [38] 3-mode tensor×in R3×3×4 in Fig. 2, where each sm×all cube
represents a scalar.
π (ζ)=(ζ γ )π (ζ) κ π (ζ),
j+1 − j j − j j−1 (7)
π 1(ζ)=0, π0(ζ)=1, j =0, ,p 1 Definition 2 (Inner Product of Two Tensors). For A,B
with − ··· − RN1×N2×···×Nd, their inner product is defined as the sum o∈f
their element-wise product
ζπ2(ζ)ρ(ζ)dζ π2 (ζ)ρ(ζ)dζ
j j+1
γj = RRRR πj2(ζ)ρ(ζ)dζ , κj+1 = RRRR πj2(ζ)ρ(ζ)dζ (8) hA,Bi=i1X,···,idA(i1,···id)B(i1,···id). (11)
Definition 3 (Frobenius Norm of A Tensor). For A
laenaddiκng0 c=oef1fi,cwiehnetr1e.πTjh(eζfi)risst pa+de1gureneiv-jaripaotelybnaosmisiaflunwcittihonas RN1×N2×···×Nd, its Frobenius norm is defined as ∈
can be obtained by normalization: A = A,A . (12)
k kF h i
p
φj(ζ)= √κ0πκj1(ζ) κj, forj =0,1,··· ,p. (9) DReNfi1n×i·t··i×onNd4is(rRaannkko-nOenifeitTceannsboersw).ritAtend-amsothdeeotuetnesroprroAduc∈t
···
of d vectors
The parameters κ ’s and γ ’s can be further used to form a
j j
symmetric tridiagonal matrix J R(p+1)×(p+1): A=v(1) v(2) v(d), withv(k) RNk (13)
∈ ◦ ···◦ ∈
γj 1, if j =k where denotes the outer product operation. This means that
− ◦
J(j,k)= √κj, if k =j+1 for1 j,k p+1. (10) d
 √κk, if k =j−1 ≤ ≤ A(i1, ,id)= v(k)(ik)forall1 ik Nk. (14)
0, otherwise ··· ≤ ≤
Let J=UΣUT be an eigenvalue decomposition, where U is Here v(k)(i ) denotekYs=t1he i -th element of vector v(k).
k k
a unitary matrix. The j-th quadrature point and weight of ζ
are Σ(j,j) and (U(1,j))2, respectively [43]. Definition 5 (Tensor Rank). The rank of A ∈ RN1×···×Nd
is the smallest positive integer r¯, such that
Challenges in High Dimension. When d is large, it is
i
difficult to implement hierarchical uncertainty quantification. r¯
First,itisnon-trivialtoobtainageneralizedpolynomialchaos A= vj(1)◦vj(2)···◦vj(d), withvj(k) ∈RNk. (15)
expansion for y , since a huge number of basis functions and Xj=1
i
It is attractive to perform tensor decomposition: given a A. ANOVA and Anchored ANOVA Decomposition
small integer r < r¯, approximate A RN1×···×Nd by a 1) ANOVA: With ANOVA decomposition [34], [64], y can
∈
rank-r tensor. Popular tensor decomposition algorithms in- be written as
clude canonical decomposition [44]–[46] and Tuker decom- y =g(ξ~)= gs(ξ~s), (18)
position [47], [48]. Canonical tensor decomposition aims to sX
approximate A by the sum of r rank-1 tensors [in the form ⊆I
where s is a subset of the full index set = 1,2, ,d .
of (15)] while minimizing the approximation error, which is Let ¯s be the complementary set of s sucIh tha{t s ·¯s··= }
normally implemented with alternating least square [45]. This and s ¯s = , and let s be the number of eleme∪nts in sI.
irdseepcirolelm-speponostseiadtiotfenonrssocdral≥ebsy3wae[sl4lm9w]a.liltThcuocthrkeeertdedinmescoeornmsaiponondsailstiiotoymneda,immbausttritixot Wξ~sh=en[ξ∩si1=,··(cid:8)·i1,∅,ξ·i|·s·|],∈i|sΩ|(cid:9)s|6=a|n∅d,hwaevesethteΩLse=beΩsgiu1e⊗m·e·a·s⊗urΩei|s|,
factors [47], [48]. This decomposition is based on singular dµ(ξ~ )= (ρ (ξ )dξ ). (19)
s¯ k k k
value decomposition. It is robust, but the number of elements
kYs¯
in the core tensor still grows exponentially with d. ∈
Then, gs(ξ~s) in ANOVA decomposition (18) is defined recur-
Alternatively, tensor-train decomposition [39]–[41] approx-
sively by the following formula
imates A RN1×···×Nd by a low-rank tensor Aˆ with
∈
E g(ξ~) = g(ξ~)dµ(ξ~)=g , if s =
0
Aˆ(i1,···id)=G1(:,i1,:)G2(:,i1,:)···Gd(:,id,:). (16) gs(ξ~s)= gˆs(cid:16)(ξ~s) (cid:17) ΩRgt(ξ~t), if s = . ∅ (20)
−t s 6 ∅
HseecroendGkind∈exRirkk,−1G×kN(:k,×ikrk,,:)anRdrkr−01×=rkrdbe=com1.esBaymfixaitnrigx t(hoer Here E is the expectatioP⊂n operator, gˆs(ξ~s) = g(ξ~)dµ(ξ~¯s),
vector when k equals 1 or∈d). To some extent, tensor-train ΩR¯s
and the integration is computed for all elements except those
decomposition have the advantages of both canonical tensor in ξ~s. From (20), we have the following intuitive results:
decomposition and Tuker decomposition: it is robust since
g is a constant term;
0
each core tensor is obtained by a well-posed low-rank matrix • if s= j , then gˆs(ξ~s) = gˆ j (ξj), gs(ξ~s) = g j (ξj) =
decomposition[39]–[41];itscaleslinearlywithdsincestoring • gˆ (ξ{)} g ; { } { }
bNalolkucn=odre(cid:15)N,tetnhasenodrtesrnrkseoq=ruitrrreasifnoordnelkyco=Om(p1No,sr·i·2t·ido),ndmi−nem(11o.6rGy) eiivfneswnueraeanssseurmroer • ibgf{soj(st}hξ~=sgˆ){jsj=(,ξ~−ksgˆ}){jaa0,knn}dd(ξgjjs,(<ξξ~ksk)),−atrhege{nsj}gˆ-(vsξ(ajξ~r)isa−)b=lge{fkgˆu}{n(j,ξckkt}i)o(ξn−js,,gξa0kn;)datnhde
• | |
decomposition (18) has 2d terms in total.
A Aˆ ε A (17)
(cid:13) − (cid:13)F ≤ k kF Example1. Considery =g(ξ~)=g(ξ ,ξ ).Since = 1,2 ,
(cid:13) (cid:13) 1 2 I { }
(cid:13) (cid:13) its subset includes , 1 , 2 and 1,2 . As a result, there
while keeping r ’s as small as possible [39].
k ∅ { } { } { }
exist four terms in the ANOVA decomposition (18):
Definition 6 (TT-Rank). In tensor-train decomposition (16) for s= , g (ξ~ )=E g(ξ~) =g is a constant;
G[rk0,r∈1,·R··r,k−rd1×] Niskc×arlkledfoTrT-kran=k. 1,···d. The vector ~r = •• forgs(ξ~)=ρ∅{(ξ1}∅),dgξ∅{1}is(ξa1)u(cid:16)=nigˆv{a1r}(cid:17)i(aξt1e)f−u0ncgt0i,onanodf ξgˆ{;1}(ξ1) =
2 2 2 1
Recently, tensor decomposition has shown promising appli- ΩR2
for s = 2 , g (ξ )=gˆ (ξ ) g , and gˆ (ξ ) =
cationsinhigh-dimensionaldataandimagecompression[50]– • g(ξ~)ρ {(ξ})dξ{2}is a2univ{a2r}iat2e f−unct0ion of ξ{;2} 2
[53], and in machine learning [54], [55]. In the uncertainty 1 1 1 2
quantification community, some efficient high-dimensional ΩfRo1r s= 1,2 , g (ξ ,ξ )=gˆ (ξ ,ξ ) g (ξ )
1,2 1 2 1,2 1 2 1 1
stochastic PDE solvers have been developed based on canon- • g (ξ{) g}. S{ince}s¯= , we h{ave}gˆ (ξ−,ξ{)}=g(ξ~−),
2 2 0 1,2 1 2
ical tensor decomposition [56]–[58] (which is called “Proper w{hi}ch is−a bi-variate fun∅ction. { }
Generalized Decomposition” in some papers) and tensor-
train decomposition [59]–[62]. In [63], a spectral tensor- Since all terms in the ANOVA decomposition are mutually
traindecompositionisproposedforhigh-dimensionalfunction orthogonal [34], [64], we have
approximation. Var g(ξ~) = Var gs(ξ~s) (21)
(cid:16) (cid:17) sX (cid:16) (cid:17)
⊆I
where Var( ) denotes the variance over the whole parameter
III. ANOVA-BASEDSURROGATEMODELEXTRACTION •
space Ω. What makes ANOVA practically useful is that for
many engineering problems, g(ξ~) is mainly influenced by the
In order to accelerate the low-level simulation, this section
terms that depend only on a small number of variables, and
developsasparsestochasticcircuit/MEMSsimulatorbasedon
thus it can be well approximated by a truncated ANOVA
anchored ANOVA (analysis of variance). Without of loss of
generality, let y =g(ξ~) denote the output of a subsystem. We decomposition
assume that y is a smooth function of the random parameters g(ξ~) gs(ξ~s), s (22)
ξ~∈Ω⊆Rd that describe the process variations. ≈|s|X≤deff ⊆I
where d d is called the effective dimension. g(ξ , ,ξ ) = g(λ (u ), ,λ (u )) = ψ(~u) with ~u =
eff 1 d 1 1 d d
(cid:28) ··· ···
[u , ,u ]. Following (25), we have
Example 2. Consider y = g(ξ~) with d = 20. In the full 1 ··· d
ANOVA decomposition (18), we need to compute over 106 p , if k ¯s
terms, which is prohibitively expensive. However, if we set ψˆs(~us)=ψ(u˜), withu˜k =(cid:26) ukk, othe∈rwise, (26)
d =2, we have the following approximation
eff where p~ = [p , ,p ] is the anchor point for ~u. The above
1 d
···
20 result can be rewritten as
g(ξ~)≈g0+Xj=1gj(ξj)+1≤jX<k≤20gj,k(ξj,ξk) (23) gˆs(ξ~s)=g(cid:16)ξ˜(cid:17),withξ˜k =(cid:26) λλkk((qξkk)),, oifthke∈rw¯sise, (27)
which contains only 221 terms.
from which we can obtain gs(ξ~s) defined in (20). Conse-
Unfortunately, it is still expensive to obtain the truncated quently, the decomposition for g(ξ~) can be obtained by using
ANOVA decomposition (22) due to two reasons. First, the ~q =[λ (p ), ,λ (p )] as an anchor point of ξ~.
1 1 d d
···
high-dimensional integrals in (20) are expensive to compute. Anchorpointselection.Itisisimportanttoselectaproper
Second, the truncated ANOVA decomposition (22) still con- anchor point [36]. In circuit and MEMS applications, we find
tains lots of terms when d is large. In the following, we that ~q =E(ξ~) is a good choice.
introduce anchored ANOVA that solves the first problem. The 2) Adaptive Implementation: In order to further reduce
second issue will be addressed in Section III-B. the computational cost, the truncated ANOVA decomposition
2) Anchored ANOVA: In order to avoid the expensive (22) can be implemented in an adaptive way. Specifically, in
multidimensional integral computation, [34] has proposed an practical computation we can ignore those terms that have
efficient algorithm which is called anchored ANOVA in [33], small variance values. Such a treatment can produce a highly
[35], [36]. Assuming that ξk’s have standard uniform distri- sparse generalized polynomial-chaos expansion.
butions, anchored ANOVA first chooses a deterministic point For a given effective dimension d d, let
eff
called anchored point ~q = [q , ,q ] [0,1]d, and then (cid:28)
replaces the Lebesgue measure1w·i·th· thde D∈irac measure k = s s , s =k , k =1, deff (28)
S { | ⊂I | | } ···
contain the initialized index sets for all k-variate terms in
dµ(ξ~ )= (δ(ξ q )dξ ). (24)
s¯ k− k k the ANOVA decomposition. Given an anchor point ~q and a
kYs¯
∈ threshold σ, starting from k=1, the main procedures of our
As a result, g =g(~q), and ANOVA-based stochastic simulator are summarized below:
0
1) Compute g , which is a deterministic evaluation;
q , if k ¯s 0
gˆs(ξ~s)=g(cid:16)ξ˜(cid:17), withξ˜k =(cid:26) ξkk, othe∈rwise. (25) 2) Fgso(rξ~esv)erbyysst∈ocShka,stciocmtepsuttiength.eTlhoew-idmimpoerntasniocnealoffugnsc(tiξ~osn)
Here ξ˜ denotes the k-th element of ξ˜ Rd, q is a fixed is measured as
k k
∈
deterministic value, and ξk is a random variable. Anchored Var gs ξ~s
ANOVAwasfurtherextendedtoGaussianrandomparameters θs = (cid:16) (cid:16) (cid:17)(cid:17) . (29)
k
in [35]. In [33], [36], [37], this algorithm was combined with Var g˜s ξ~˜s
stochastic collocation to efficiently solve high-dimensional jP=1˜sP∈Sj (cid:16) (cid:16) (cid:17)(cid:17)
stochastic partial differential equations.
3) Updatetheindexsetsifθs <σ fors k.Specifically,
~gqEˆx2=am(ξ[qp21)l,eq=23].,gCw(qoe1n,hsξiad2ve)eragyn0=d=ggˆ(ξg11(,,2qξ12(,)ξq.12,W)ξ,2itg)ˆh{1=a}n(ξg1a()nξc1=h,oξ2rge)(d.ξ1Cp,ooq2imn)-t, cOfoonrnctkeais<n0sijasl≤lreemdleeomffv,eednw,tsewocefhdesoc,ktnhoiettsnniwenedederxteo∈mseeSovtvaselu0sa∈0tefSrgojs.m0(Iξ~fSss0j)0.
p{uti}ng these quantities does n{ot i}nvolve any high-dimensional in the subsequent computation.
integrations. 4) Set k=k+1, and repeat steps 2) and 3) until k =d .
def
Example 4. Let y=g(ξ~), ξ~ R20 and d = 2. Anchored
eff
∈
B. Adaptive Anchored ANOVA for Circuit/MEMS Problems ANOVA starts with
1) ExtensiontoGeneralCases: InmanycircuitandMEMS = j and = j,k .
S1 {{ }}j=1, ,20 S2 {{ }}1 j<k 20
problems, the process variations can be non-uniform and non- ··· ≤ ≤
Gaussian. We show that anchored ANOVA can be applied to For k=1, we first utilize stochastic testing to calculate gs(ξ~s)
such general cases. and θs for every s 1. Assume
∈S
Observation: The anchored ANOVA in [34] can be applied
θ >σ, θ >σ, andθ <σ forallj >2,
1 2 j
if ρk(ξk)>0 for any ξk Ωk. { } { } { }
∈
Proof: Let u denote the cumulative density function for implying that only the first two parameters are important to
k
ξk, then uk can be treated as a random variable uniformly the output. Then, we only consider the coupling of ξ1 and ξ2
distributed on [0,1]. Since ρk(ξk) > 0 for any ξk ∈ Ωk, in S2, leading to
there exists ξ = λ (u ) which maps u to ξ . Therefore, = 1,2 .
k k k k k 2
S {{ }}
Algorithm 1 Stochastic Testing Circuit/MEMS Simulator 3) Global Sensitivity Analysis: Since each term gs(ss) is
Based on Adaptive Anchored ANOVA. computedbystochastictesting,Algorithm1providesasparse
1: Initialize k’s and set β =0; generalized polynomial-chaos expansion for the output of
S
2: At the anchor point, run a deterministic circuit/MEMS interest: y= y H (ξ~), where most coefficients are zero.
α~ α~
simulation to obtain g0, and set y =g0; From this re|α~sPu|≤ltp, we can identify how much each parameter
3: for k =1, , deff do
4: for each·s·· k do contributes to the output by global sensitivity analysis. Two
∈S kinds of sensitivity information can be used to measure the
5: runstochastictestingsimulatortogetthegeneralized
polynomial-chaos expansion of gˆs(ξ~s) ; importance of parameter ξk: the main sensitivity Sk and total
sensitivity T , as computed below:
6: get the generalized polynomial-chaos expansion of k
7: gusp(dξ~ast)eaβcc=orβdi+ngVtoar(2(cid:16)0g)s;(ξ~s)(cid:17); Sk = αk6=0,PVαja6=rk(=y0)|yα~|2, Tk = αkPV6=a0r|(yyα~)|2. (31)
8: update y =y+gs(ξ~s);
9: end for IV. ENABLINGHIGH-LEVELSIMULATIONBY
10: for each s k do TENSOR-TRAINDECOMPOSITION
∈S
11: θs =Var gs(ξ~s) /β; In this section, we show how to accelerate the high-level
(cid:16) (cid:17)
12: if θs <σ non-Monte-Carlo simulation by handling the obtained high-
13: for any index set s0 j with j >k, remove dimensional surrogate models with tensor-train decomposi-
s0 from j if s s0.∈S tion [39]–[41].
S ⊂
14: end if
15: end for A. Tensor-Based Three-Term Recurrence Relation
16: end for
In order to obtain the orthonormal polynomials and Gauss
quadrature points/weights of ζ, we must implement the three-
term recurrence relation in (7). The main bottleneck is to
Consequently, for k = 2 we only need to calculate one bi-
compute the integrals in (8), since the probability density
variate function g (ξ ,ξ ), yielding
{1,2} 1 2 function of ζ is unknown.
g ξ~ g + g ξ~ + g ξ~ For simplicity, we rewrite the integrals in (8) as E(q(ζ)),
0 s s s s
(cid:16) (cid:17)≈ 2s0P∈S1 (cid:16) (cid:17) sP∈S2 (cid:16) (cid:17) wdeinthsitqy(ζfu)n=ctioφn2j(oζf)ζorisqn(oζt)g=iveζnφ,2jw(ζe)c.oSminpcueteththeepirnotebgarbaillitiyn
=g + g (ξ )+g (ξ ,ξ ).
0 jP=1 {j} j {1,2} 1 2 the parameter space Ω:
Thepseudocodesofourimplementationaresummarizedin E(q(ζ))= q f ξ~ ρ(ξ~)dξ dξ , (32)
1 d
Alg.1.Lines10to15showshowtoadaptivelyselecttheindex Z (cid:16) (cid:16) (cid:17)(cid:17) ···
sets. Let the final size of be and the total polynomial Ω
k k
order in the stochastic teSsting s|iSmu|lator be p, then the total where f(ξ~) is a sparse generalized polynomial-chaos expan-
number of samples used in Alg. 1 is sion for ζ obtained by
deff (k+p)! ζ =f(ξ~)= (y−E(y)) = yˆ H (ξ~). (33)
N =1+kX=1|Sk| k!p! . (30) pVar(y) |α~X|≤p α~ α~
NotethatallunivariatetermsinANOVA(i.e., s =1)arekept We compute the integral in (32) with the following steps:
inourimplementation.FormostcircuitandM|E|MSproblems, 1) We utilize a multi-dimensional Gauss quadrature rule:
setting the effective dimension as 2 or 3 can achieve a high m1 md d
accuracy due to the weak couplings among different random E(q(ζ)) q f ξi1, ,ξid wik
≈ ··· 1 ··· d k
parameters.Formanycases,theunivariatetermsdominatethe iX1=1 iXd=1 (cid:0) (cid:0) (cid:1)(cid:1)kY=1
output of interest, leading to a near-linear complexity with (34)
respect to the parameter dimensionality d. where mk is the number of quadrature points for ξk,
Remarks. Anchored ANOVA works very well for a large (ξkik,wkik) denotes the ik-th Gauss quadrature point and
class of MEMS and circuit problems. However, in practice weight.
we also find a small number of examples (e.g., CMOS ring 2) Wedefinetwod-modetensorsQ,W Rm1×m2···×md,
∈
oscillators) that cannot be solved efficiently by the proposed with each element defined as
algorithm, since many random variables affect significantly Q(i , i )=q f ξi1, ,ξid ,
1 ··· d 1 ··· d
the output of interest. For such problems, it is possible (cid:0)d (cid:0) (cid:1)(cid:1) (35)
W(i , i )= wik,
to reduce the number of dominant random variables by a 1 ··· d k
lineartransform[65]beforeapplyinganchoredANOVA.Other kQ=1
for 1 i m . Now we can rewrite (34) as the inner
techniques such as compressed sensing can also be utilized to ≤ k ≤ k
product of Q and W:
extract highly sparse surrogate models [66]–[69] in the low-
level simulation of our proposed hierarchical framework. E(q(ζ)) Q,W . (36)
≈h i
For simplicity, we set m =m in this manuscript. Fast evaluation of Q(i , ,i ). In order to reduce the
k 1 d
• ···
The cost of computing the tensors and the tensor inner costofevaluatingQ(i1, ,id),wefirstconstructalow-
productisO(md),whichbecomesintractablewhendislarge. rank tensor train Aˆ for··t·he intermediate-level random
Fortunately, both Q and W have low tensor ranks in our parameter ζ, such that
applications, and thus the high-dimensional integration (32)
A Aˆ ε A , A(i , ,i )=f ξi1, ,ξid .
can be computed very efficiently in the following way: (cid:13) − (cid:13)F ≤ k kF 1 ··· d 1 ··· d
1) Low-rank representation of W. W can be written as O(cid:13)(cid:13)nce Aˆ(cid:13)(cid:13)is obtained, Q(i , ,i ) can be e(cid:0)valuated by (cid:1)
1 d
a rank-1 tensor ···
W =w(1)◦w(2)···◦w(d), (37) Q(i1,··· ,id)≈q(cid:16)Aˆ(i1,··· ,id)(cid:17), (41)
where w(k) = [w1; ;wm] Rm 1 contains all whichreducestoacheaplow-orderunivariatepolynomial
k ··· k ∈ × evaluation. However, computing Aˆ(i , ,i ) by di-
Gaussquadratureweightsforparameterξ .Clearly,now 1 d
k ···
we only need O(md) memory to store W. rectlyevaluatingA(i1, ,id)inTT crosscanbetime-
2) Low-rank approximation for Q. Q can be well ap- consuming, since ζ =·f·(·ξ~) involves many multivariate
proximated by Qˆ with high accuracy in a tensor-train basis functions.
Fast evaluation of A(i , ,i ). The evaluation of
format [39]–[41]: 1 d
• ···
A(i , ,i ) can also be accelerated by exploiting the
1 d
Qˆ(i1, id)=G1(:,i1,:)G2(:,i1,:) Gd(:,id,:) special·s·t·ructureoff(ξ~).Itisknownthatthegeneralized
··· ···
(38) polynomial-chaos basis of ξ~ is
with a pre-selected error bound (cid:15) such that
d
(cid:13)Q−Qˆ(cid:13)F ≤εkQkF . (39) Hα~ (cid:16)ξ~(cid:17)=kY=1ϕα(kk)(ξk), α~ =[α1,··· ,αd] (42)
(cid:13) (cid:13)
For many circu(cid:13)it and M(cid:13)EMS problems, a tensor train
whereϕ(k)(ξ )isthedegree-α orthonormalpolynomial
withverysmallTT-rankscanbeobtainedevenwhen(cid:15)= αk k k
of ξ , with 0 α p. We first construct a 3-mode
10−12 (which is very close to the machine precision). tensokr X Rd≤(p+k1)≤m indexed by (k,α +1,i ) with
3) Fast computation of (36). With the above low-rank ∈ × × k k
tensor representations, the inner product in (36) can be X (k,α +1,i )=ϕ(k) ξik (43)
accurately estimated as k k αk k
(cid:0) (cid:1)
where ξik is the i -th Gauss quadrature point for pa-
m k k
Qˆ,W =T T , withT = wikG (:,i ,:) rameter ξk [as also used in (34)]. Then, each element of
D E 1··· d k iXk=1 k k k A(i1,··· ,id) can be calculated efficiently as
(40)
d
Now the cost of computing the involved high-
A(i , ,i )= ~y X (k,α +1,i ) (44)
1 d α~ k k
dimensional integration dramatically reduces to ···
α~X<p kY=1
O(dmr2), which only linearly depends the parameter | |
dimensionality d. without evaluating the multivariate polynomials. Con-
structing X does not necessarily need d(p+1)m poly-
nomial evaluations, since the matrix X (k,:,:) can be
B. Efficient Tensor-Train Computation
reused for any other parameter ξ that has the same type
j
Now we discuss how to obtain a low-rank tensor train. of distribution with ξ .
k
An efficient implementation called TT cross is described
In summary, we compute a tensor-train decomposition for
in [41] and included in the public-domain MATALB package
Q as follows: 1) we construct the 3-mode tensor X defined
TT Toolbox [70]. In TT cross, Skeleton decomposition is in (43); 2) we call TT cross to compute Aˆ as a tensor-
utilized to compress the TT-rank rk by iteratively searching a train decomposition of A, where (44) is used for fast element
rank-rk maximum-volume submatrix when computing Gk. A evaluation; 3) we call TT cross again to compute Qˆ, where
major advantage of TT cross is that we do not need to know
(41) is used for the fast element evaluation of Q. With the
Q a-priori. Instead, we only need to specify how to evaluate
abovefasttensorelementevaluations,thecomputationtimeof
the element Q(i , ,i ) for a given index (i , ,i ). As
1 d 1 d TT cross can be reduced from dozens of minutes to several
··· ···
shown in [41], with Skeleton decompositions a tensor-train
seconds to generate some accurate low-rank tensor trains for
decomposition needs O(ldmr2) element evaluations, where l
our high-dimensional surrogate models.
is the number of iterations in a Skeleton decomposition. For
example, when l = 10, d = 50, m = 10 and r = 4 we may
need up to 105 element evaluations, which can take about C. Algorithm Summary
one hour since each element of Q is a high-order polynomial Given the Gauss quadrature rule for each bottom-level
function of many bottom-level random variables ξ~. random parameter ξ , our tensor-based three-term recurrence
k
In order to make the tensor-train decomposition of Q fast, relation for an intermediate-level random parameter ζ is sum-
we employ some tricks to evaluate more efficiently each marized in Alg. 2. This procedure can be repeated for all ζ ’s
i
element of Q. The details are given below. to obtain their univariate generalized polynomial-chaos basis
Algorithm2Tensor-basedgeneralizedpolynomial-chaosbasis Vdd
and Gauss quadrature rule construction for ζ.
1: Initialize: φ0(ζ) = π0(ζ) = 1, φ1(ζ) = π1(ζ) = ζ, κ0 = L R L Parameter Value
κ1 =1, γ0 =0, a=1; Vdd 2.5V
2: Compute a low-rank tensor train Aˆ for ζ; C
3: oCbotmaipnuγte a=lowQˆ-,raWnk tevniaso(r40tr)a;in Qˆ for q(ζ) = ζ3, and V1 Cm Vctrl Cm V2 VCctrl 100-2.8.53VfF
4: for j =12, D, p doE Cm Cm R 110Ω
5: get πj(ζ·)·=· (ζ γj 1)πj 1(ζ) κj 1πj 2(ζ) ; L 8.8nH
6: construct a low-−rank−tenso−r train−Qˆ f−or q(−ζ)=πj2(ζ), (W/L)n 4/0.25
and compute aˆ= Qˆ,W via (40) ; Mn Mn Iss 0.5mA
7: κj =aˆ/a, and updDate a=Eaˆ ; Rss 106Ω
8: constructalow-ranktensortrainQˆ forq(ζ)=ζπ2(ζ),
j R I
ss ss
and compute γ = Qˆ,W /a ;
j
D E
190:: endnoformralization: φj(ζ)= √πκj0(·ζ··)κj ; FasigC.3m.),Swcihthem1a8t4icroafndthoemospcairlalamtoertecrisrciunittowtaitlh. 4MEMScapacitors(denoted
11: Form matrix J in (10);
12: Eigenvalue decomposition: J=UΣUT ;
13: Compute the Gauss-quadrature abscissa ζj =Σ(j,j) and RF Conducting
weight wj =(U(1,j))2 for j =1, ,p+1 ; Path
···
functions and Gauss quadrature rules, and then the stochastic
testing simulator [14]–[16] (and any other standard stochastic Secondary RF Capacitor
spectral method [7]–[9]) can be employed to perform high- Actuator Contact Bump Primary
Actuator
level stochastic simulation.
Remarks. 1) If the outputs of a group of subsystems are
identically independent, we only need to run Alg. 2 once
Fig.4. 3-DschematicoftheRFMEMScapacitor.
and reuse the results for the other subsystems in the group.
2) When there exist many subsystems, our ANOVA-based
stochastic solver may also be utilized to accelerate the high- shown in [73], the performance of this MEMS switch can be
level simulation. influenced significantly by process variations.
In our numerical experiments, we use 46 independent ran-
dom parameters with Gaussian and Gamma distributions to
V. NUMERICALRESULTS
describe the material (e.g, conductivity and dielectric con-
In this section, we verify the proposed algorithm on a stants), geometric (e.g., thickness of each layer, width and
MEMS/IC co-design example with high-dimensional random length of each mechanical component) and environmental
parameters.AllsimulationsareruninMATLABandexecuted (e.g., temperature) uncertainties of each switch. For each
on a 2.4GHz laptop with 4GB memory. random parameter, we assume that its standard deviation is
3% of its mean value. In the whole circuit, we have 184
random parameters in total. Due to such high dimensionality,
A. MEMS/IC Example
simulating this circuit by stochastic spectral methods is a
In order to demonstrate the application of our hierarchical challenging task.
uncertainty quantification in high-dimensional problems, we In the following experiments, we simulate this challenging
consider the oscillator circuit shown in Fig. 3. This oscillator designcaseusingourproposedhierarchicalstochasticspectral
has four identical RF MEMS switches acting as tunable ca- methods.Wealsocompareouralgorithmwithothertwokinds
pacitors.TheMEMSdeviceusedinthispaperisaprototyping of hierarchical approaches listed in Table I. In Method 1,
model of the RF MEMS capacitor reported in [71], [72]. bothlow-levelandhigh-level simulationsuseMonteCarlo,as
Since the MEMS switch has a symmetric structure, we suggested by [32]. In Method 2, the low-level simulation uses
construct a model for only half of the design, as shown in our ANOVA-based sparse simulator (Alg. 1), and the high-
Fig. 4. The simulation and measurement results in [42] show level simulation uses Monte Carlo.
that the pull-in voltage of this MEMS switch is about 37 V.
When the control voltage is far below the pull-in voltage,
B. Surrogate Model Extraction
the MEMS capacitance is small and almost constant. In this
paper,wesetthecontrolvoltageto2.5V,andthustheMEMS In order to extract an accurate surrogate model for the
switch can be regarded as a small linear capacitor. As already MEMS capacitor, Alg. 1 is implemented in the commercial
TABLEII
SURROGATEMODELEXTRACTIONWITHDIFFERENTσVALUES.
σ #|s|=1 #|s|=2 #|s|=3 #ANOVAterms #nonzerogPCterms #samples
0.5 46 0 0 47 81 185
0.1to10−3 46 3 0 50 90 215
10−4 46 10 1 58 112 305
10−5 46 21 1 69 144 415
TABLEI
DIFFERENTHIERARCHICALSIMULATIONMETHODS. 0.4
main sensitivity
total sensitivity
Method Low-levelsimulation High-levelsimulation
0.3
Proposed Alg.1 stochastictesting[15]
Method1[32] MonteCarlo MonteCarlo
Method2 Alg.1 MonteCarlo 0.2
1.4 0.1
surrogate (s =10−2)
1.2 Monte Carlo (5000)
0
0 5 10 15 20 25 30 35 40 45
1 Parameter index
0.8
Fig. 6. Main and total sensitivities of different random parameters for the
0.6
RFMEMScapacitor.
0.4
0.2 9
TT−rank
8
0
5.5 6 6.5 7 7.5 8 8.5 9 7
MEMS capacitor (fF)
6
5
Fig.5. Comparisonofthedensityfunctionsobtainedbyoursurrogatemodel
andby5000-sampleMonteCarloanalysisoftheoriginalMEMSequation. 4
3
2
1
network-basedMEMSsimulationtoolMEMS+[74]ofCoven-
tor Inc. Each MEMS switch is described by a stochastic 00 5 10 15 20 25 30 35 40 45 50
Parameter index
differential equation [c.f. (1)] with consideration of process
variations. In order to compute the MEMS capacitor, we can Fig.7. TT-rankforthesurrogatemodeloftheRFMEMScapacitor.
ignore the derivative terms and solve for the static solutions.
By setting σ =10 2, our ANOVA-based stochastic MEMS
−
changing the threshold σ. Table II has listed the number of
simulatorgeneratesasparse3rd-ordergeneralizedpolynomial
obtained ANOVA terms, the number of non-zero generalized
chaos expansion with only 90 non-zero coefficients, requiring
polynomial chaos (gPC) terms and the number of required
only 215 simulation samples and 8.5 minutes of CPU time
simulation samples for different values of σ. From this table,
in total. This result has only 3 bivariate terms and no three-
we have the following observations:
variabletermsinANOVAdecomposition,duetotheveryweak
couplings among different random parameters. Setting σ = 1) Whenσ islarge,only46univariateterms(i.e.,theterms
10 2 can provide a highly accurate generalized polynomial with s =1) are obtained. This is because the variance
− | |
chaosexpansionfortheMEMScapacitor,whichhasarelative of all univariate terms are regarded as small, and thus
erroraround10 6 (intheL sense)comparedtothatobtained all multivariate terms are ignored.
− 2
by setting σ =10 5. 2) Whenσ isreduced(forexample,to0.1),threedominant
−
By evaluating the surrogate model and the original model bivariate terms (with s = 2) are included by consid-
| |
(by simulating the original MEMS equation) with 5000 sam- ering the coupling effects of the three most influential
ples, we have obtained the same probability density curves random parameters. Since the contributions of other
shown in Fig. 5. Note that using the standard stochastic parameters are insignificant, the result does not change
testingsimulator[14]–[16]requires18424basisfunctionsand even if σ is further decreased to 10−3.
simulation samples for this high-dimensional example, which 3) A three-variable term (with s =3) and some bivariate
| |
is prohibitively expensive on a regular computer. When the coupling terms among other parameters can only be
effective dimension deff is set as 3, there should be 16262 captured when σ is reduced to 10−4 or below. In this
terms in the truncated ANOVA decomposition (22). However, case, the effect of some non-dominant parameters can
duetotheweakcouplingsamongdifferentrandomparameters, be captured.
only 90 of them are non-zero. Fig. 6 shows the global sensitivity of this MEMS capacitor
We can get surrogate models with different accuracies by with respect to all 46 random parameters. The output is
(a) (b)
0.8 10
s
0.6 n 5
o
cti
s n
weight0.4 basis fu 0 k=0
C k=1
0.2 P −5
g k=2
k=3
0 −10
−4 −2 0 2 4 −4 −2 0 2 4
Gauss quadrature points z
Fig.8. (a)Gaussquadratureruleand(b)generalizedpolynomialchaos(gPC)basisfunctionsfortheRFMEMScapacitor.
dominated by only 3 parameters. The other 43 parameters 4) Map ~z(τ,ζ~) to the original time axis, we obtain the
contributetoonly2%ofthecapacitor’svariance,andthustheir periodic steady state of ~x(t,ζ).
main and total sensitivities are almost invisible in Fig. 6. This In order to apply stochastic testing in Step 3), we need to
explains why the generalized polynomial-chaos expansion is computesomespecializedorthonormalpolynomialsandGauss
highly sparse. Similar results have already been observed in quadraturepointsforeachintermediate-levelparameterζ .We
i
the statistical analysis of CMOS analog circuits [1]. use 9 quadrature points for each bottom-level parameter ξ to
k
evaluate the high-dimensional integrals involved in the three-
termrecurrencerelation.Thisleadsto946 functionevaluations
C. High-Level Simulation
at all quadrature points, which is prohibitively expensive.
The surrogate model obtained with σ = 10 2 is imported To handle the high-dimensional MEMS surrogate models,
−
into the stochastic testing circuit simulator described in [14]– the following tensor-based procedures are employed:
[16] for high-level simulation. At the high-level, we have the With Alg. 2, a low-rank tensor train of ζ is first con-
1
•
following differential equation to describe the oscillator: structed for an MEMS capacitor. For most dimensions
the rank is only 2, and the highest rank is 4, as shown in
d~q ~x(t,ζ~),ξ~ Fig. 7.
(cid:16) (cid:17) +f~ ~x(t,ζ~),ζ~,u =0 (45)
Using the obtained tensor train, the Gauss quadrature
dt
(cid:16) (cid:17) •
points and generalized polynomial chaos basis functions
where the input signal u is constant, ζ~=[ζ1, ,ζ4] R4 are efficiently computed, as plotted in Fig. 8.
··· ∈
are the intermediate-level random parameters describing the
The total CPU time for constructing the tensor trains
four MEMS capacitors. Since the oscillation period T(ζ~) now and computing the basis functions and Gauss quadrature
depends on the MEMS capacitors, the periodic steady-state
points/weights is about 40 seconds in MATALB. If we di-
can be written as ~x(t,ζ) = ~x(t+T(ζ~),ζ). We simulate the rectly evaluate the high-dimensional multivariate generalized
stochastic oscillator by the following steps [15]:
polynomial-chaos expansion, the three-term recurrence rela-
1) Choose a constant T >0 to define an unknown scaling tion requires almost 1 hour. The obtained results can be
0
factor a(ζ~) = T(ζ~)/T and a scaled time axis τ = reused for all MEMS capacitors since they are independently
0
t/α(ζ~). With this scaling factor, we obtain a reshaped identical.
waveform ~z(τ,ζ~) = ~x(t/a(ζ~),ζ~). At the steady state, With the obtained basis functions and Gauss quadrature
we have ~z(τ,ζ~) = ~z(τ + T ,ζ~). In other words, the points/weights for each MEMS capacitor, the stochastic pe-
0
reshaped waveform has a period T independent of ζ~. riodic steady-state solver [15] is called at the high level to
0
2) Rewrite (45) on the scaled time axis: simulate the oscillator. Since there are 4 intermediate-level
parameters ζ ’s, only 35 basis functions and testing samples
i
d~q ~z(τ,ζ~),ξ~ are required for a 3rd-order generalized polynomial-chaos
(cid:16) (cid:17) +a(ζ~)f~ ~z(τ,ζ~),ζ~,u =0. (46) expansion, leading to a simulation cost of only 56 seconds
dτ
(cid:16) (cid:17)
in MATLAB.
3) Approximate ~z(τ,ζ~) and a(ζ~) by generalized polyno- Fig. 9 shows the waveforms from our algorithm at the
mial chaos expansions of ζ~. Then, convert (46) to a scaled time axis τ = t/a(ζ~). The high-level simulation gen-
larger-scale deterministic equation by stochastic testing. erates a generalized polynomial-chaos expansion for all nodal
Solve the resulting deterministic equation by shooting voltages, branch currents and the exact parameter-dependent
Newton with a phase constraints, which would provide period. Evaluating the resulting generalized polynomial-chaos
the coefficients in the generalized polynomial-chaos expansion with 5000 samples, we have obtained the density
expansions of ~z(τ,ζ~) and α(ζ~) [15]. functionofthefrequency, which isconsistentwiththosefrom