arXiv:1610.03863v4 [cs.CE] 9 Dec 2016

High Dimensional Uncertainty Quantification for
an Electrothermal Field Problem using Stochastic
Collocation on Sparse Grids and Tensor Train
Decompositions
D. Loukrezis1,2 , U. Römer1,2 , T. Casper1,2 , S. Schöps1,2 , and
H. De Gersem1,2
1

Technische Universität Darmstadt, Institut für Theorie Elektromagnetischer
Felder, Schloßgartenstraße 8, D-64289, Darmstadt, Germany
2
Technische Universität Darmstadt, Graduate School for Computational
Engineering, Dolivostraße 15, D-64293, Darmstadt, Germany
Abstract
The temperature developed in bondwires of integrated circuits (ICs)
is a possible source of malfunction, and has to be taken into account
during the design phase of an IC. Due to manufacturing tolerances, a
bondwire’s geometrical characteristics are uncertain parameters, and as
such their impact has to be examined with the use of uncertainty quantification (UQ) methods. Considering a stochastic electrothermal problem
featuring twelve (12) bondwire-related uncertainties, we want to quantify
the impact of the uncertain inputs onto the temperature developed during
the duty cycle of an IC. For this reason we apply the stochastic collocation (SC) method on sparse grids (SGs), which is considered the current
state-of-the-art. We also implement an approach based on the recently
introduced low-rank tensor decompositions, in particular the tensor train
(TT) decomposition, which in theory promises to break the curse of dimensionality. A comparison of both methods is presented, with respect
to accuracy and computational effort.
keywords— electrothermal field simulation, IC packaging, low-rank
tensor decompositions, sparse grids, stochastic collocation, tensor trains,
uncertainty quantification.

1

Introduction

A thorough electrothermal analysis is a key step in the design of integrated
circuits (ICs). One important aspect is the Joule heating of the bondwires
which connect the IC to its surrounding packaging, as wire fusing is a common
source of failure. Thermal heating is affected by uncertainties introduced in the
manufacturing process, e.g. with respect to the bondwire’s length. More often
1

than ever before uncertainty quantification (UQ) methods are applied, with the
aim of assessing the risk of fusing or failure. In this way, one may avoid costly
and time-consuming over-engineering.
In a stochastic setting, state-of-the-art UQ methods are the stochastic Galerkin and collocation method [1, 26]. These methods use a spectral approximation
based on a global multivariate polynomial basis. Given a high regularity of the
solution with respect to the random input variables, fast convergence rates can
be achieved. Hence, the computational cost of using stochastic Galerkin and
collocation can be affordable even for applications with a high complexity.
This might no longer be true as the number of random inputs increases.
Referring to the bondwire case, the number of wires can be large and, at the
same time, randomness hampers the exploitation of any symmetry properties,
which could otherwise be used in order to reduce the computational effort. The
difficulty arising when approximating high-dimensional problems is called the
curse of dimensionality [5]. In the case of standard Galerkin and collocation
methods, this cannot simply be overcome with a further increase of computational power, as the complexity scales exponentially with the number of random
inputs. Hence, the construction of more dedicated schemes is in order.
The efficiency of stochastic spectral methods can be enhanced for highdimensional problems using sparse approximations. In this case, dimension independent approximation rates can be derived under very strict assumptions on
the solution regularity and the input randomness [10]. Non-intrusive collocation
based on sparse grids [3, 6, 18] outperforms Monte Carlo and polynomial tensor
approximations when the underlying functions possess high (mixed) regularity.
However, if some of the assumptions are not fully satisfied, the convergence of
sparse approximations is generally not free of dimensionality problems. A recent application of sparse approximations in electromagnetics is described in [7]
based on the least angle regression method.
Recently, based on the observation that multivariate polynomial bases have
an underlying tensor structure, low-rank tensor approximations have gained increasing interest [13, 14, 15, 16]. First results on low-rank tensor approximations
show an improved asymptotic convergence rate with respect to the number of
random inputs compared to sparse grids [24], yet many theoretical questions
persist. Moreover, the computational effort for more complicated examples has
hardly been investigated. It is the aim of the present paper to contribute in this
direction.
This paper compares the accuracy and computational complexity of stochastic
collocation on sparse grids and the one of tensor train decomposition. Both
methods are considered in the context of electrothermal field simulations, e.g.
occurring in IC packaging problems, characterized by high dimensional uncertainty, i.e. with a comparably large number of uncertain parameters.
The rest of the paper is structured as follows: Section 2 deals with an electrothermal field problem and its numerical discretization. In Section 3 the considered uncertainties and quantities of interest are described. Section 4 describes
both stochastic collocation on sparse grids and low-rank tensor approximations
using tensor trains. A comparison of both methods when applied to the con2

sidered electrothermal problem is presented in Section 5.

2

Electro-thermal Field Problem

Due to the Joule heating effect, the temperature of a wire increases if an electrical current is applied to it. In turn, temperature changes affect the wire’s
material parameters. Even if the dependence of the volumetric heat capacity on
the temperature is neglected, the nonlinearity in electrical and thermal conductivities with respect to temperature remains. This electrothermal coupling effect
is presented in the following subsections, for the continuous and the discrete
case, respectively.

2.1

Continuous Problem

The electrical part of the problem can be described with the current continuity
equation. If capacitive effects are ignored, it suffices to write the electro-kinetic
problem with suitable boundary conditions
∇ · σ(·, T )∇φ = 0

in D × I,
on ΓDir,i × I,

φi = V i
σ(·, T )∇φ · ~n = 0

on

(1)

(Γ \ ΓDir ) × I,

where D denotes the computational domain, I the time interval, Γ the boundary
of the computational domain, ΓDir the Dirichlet boundary, Vi the potential
applied on the i-th contact pad, ~n the outer normal vector, σ the electrical
conductivity, φ the electric potential and T the temperature.
Thermal heat flow can be attributed to conduction, convection and radiation. The thermodynamic power balance is described with the transient heat
equation, which reads
ρc ∂t T − ∇ · λ(·, T )∇T = q(·, T, φ)

in D × I,
4

λ(·, T )∇T · ~n = hT + εσSB T − hT∞ −

4
εσSB T∞

on Γ × I,

(2)

on Γ × {0},

T = T∞ ,

where ρc denotes the volumetric heat capacity, λ the thermal conductivity, q
the power density of the system’s heat sources, h the heat transfer coefficient,
ε the emissivity, σSB the Stefan-Boltzmann constant and T∞ the ambient temperature. The heat sources comprise of the electrical contribution qel due to
Joule heating and the bondwire contribution qbw , such that q = qel + qbw . Heat
exchange with the environment is modeled as convective and radiative boundary
conditions, see [9].
The combination of the transient heat equation with the electro-kinetic problem yields a nonlinear electrothermal system, where the electrothermal coupling
is established due to the dependencies of the material parameters λ and σ on

3

the temperature T and of the heat source q on the electric potential φ. The
Joule heating due to the electro-kinetic problem reads
qel = (∇φ)> σ(T )∇φ.

2.2

(3)

Discrete Problem

The coupled electrothermal problem is discretized in space with the finite integration technique (FIT) [25] on a primal and a dual Cartesian mesh. The
discrete unknowns, i.e. the electric potentials Φ and the temperatures T, are
allocated at the nodes of the primary grid. Omitting specific implementation
details and initial conditions, the discrete analogue to (1), (2) reads
−S̃Mσ (T)GΦ = 0,

(4)

Mρc Ṫ − S̃Mλ (T)GT = Q(T, Φ),

(5)

where S̃> = −G is a topological operator, with G ∈ {−1, 0, 1}n×n being the
edge-node incidence matrix, Mσ the electrical conductance mass matrix, Mλ the
thermal conductance mass matrix, Mρc the thermal capacitance mass matrix
and Q the discrete counterpart of the heat sources.
Bondwires are not resolved by the FIT grid, but are modeled with a lumped
element approach, where pairs of mesh points are connected by an electrothermal lumped element. This approach is illustrated in Figure 1. Both conductances depend on the wire’s material properties and geometrical dimensions,
i.e.,
Gel
bw,j =

σbw,j Abw,j
`bw,j

and

Gth
bw,j =

λbw,j Abw,j
,
`bw,j

(6)

where `bw,j denotes the length of the j-th bondwire and Abw,j its cross-section
th
surface. The conductance terms Gel
bw,j , Gbw,j have to be assembled into the
system matrices of (4), (5), based on the bondwires’ connection positions. Then,
Gel
bw (T )

→
Figure 1: Bondwire lumped element modeling of the electrical conductance.

4

the extended electrothermal system reads
−S̃Mσ (T)GΦ +

N
bw
X

>
Pj Gel
bw,j (Tbw,j )Pj Φ = 0,

(7)

>
Pj Gth
bw,j (Tbw,j )Pj T = Q(T, Φ),

(8)

j=1

Mρc Ṫ − S̃Mλ (T)GT +

N
bw
X
j=1

where Nbw denotes the amount of bondwires, Pj ∈ {−1, 0, 1}n×n is a bondwire
gradient vector which handles the incidences between the contacts of the j-th
bondwire and the grid, and Tbw,j = X>
j T is the average temperature of the jth bondwire, with vector Xj being used to average the temperature from both
connection nodes of the bondwire. The bondwires’ contribution to the heat
sources reads
Qbw =

N
bw
X

Xj Qbw,j ,

(9)

j=1
>
where Qbw,j = Φ> Pj Gel
bw,j (Tbw,j )Pj Φ.

3
3.1

Uncertainties and Quantities of Interest
Input Uncertainties

The bondwire’s temperature depends on the wire’s material properties, its geometry and on the applied electric voltage. Assuming that the electric voltage
can be accurately controlled and that material properties are known with high
precision, the only remaining source of uncertainty is related to the wire’s geometry. Further assuming that a wire can be fabricated with an almost exact
cross-section, the wire’s length remains as the only uncertain parameter.
A wire’s length depends on the bonding process during the IC fabrication.
The distance between the contact pad and the chip determines the minimum
length of a wire, `min . Any deviation from the ideal position on the contact pad
adds an elongation ∆s to the minimum length of the wire. Any bending of the

`min + ∆s
∆h
`min
(a) Elongation due to misplaced contact.

(b) Elongation due to wire bending.

Figure 2: Bondwire length variations due to fabrication tolerances.

5

wire results in an additional elongation ∆h, yielding the total wire length
`tot = `min + ∆s + ∆h.

(10)

These elongations are illustrated in Figure 2. Due to the uncertainties concerning pad positions and wire bending, the wire’s length can be considered as an
uncertain parameter. Instead of the total wire length, we choose the uncertain
parameter to be the relative elongation, which is defined as
δ = (`tot − `min ) /`tot .

(11)

Specifically, each bondwire’s relative elongation is drawn from the range [0.122, 0.218]
[9] in a uniform manner. Hence, we consider one random variable (RV) Yn
per bondwire in the IC, corresponding to the random relative elongation and
following a uniform distribution, i.e. Yn ∼ U[0.122, 0.218]. RVs Yn are assumed to be stochastically independent. Then the input RVs are given as
Y = (Y1 , Y2 , ..., YNbw ), with joint probability density function (PDF) % (y).
Here, y = Y(ω) denotes a realization of the input RVs, with ω ∈ Ω being
a random outcome, such that Y : Ω → RNbw .

3.2

Output Uncertainties and Quantities of Interest

Uncertainty in the inputs of the problem also leads to uncertainty in the discrete
unknowns through the stochastic problem
−S̃Mσ (T)GΦ +

N
bw
X

>
Pj Gel
bw,j (Y, Tbw,j )Pj Φ = 0,

(12)

>
Pj Gth
bw,j (Y, Tbw,j )Pj T = Q(Y, T, Φ).

(13)

j=1

Mρc Ṫ − S̃Mλ (T)GT +

N
bw
X
j=1

In turn, the uncertainty in the solution has to be taken into consideration for
the estimation of quantities of interest (QoI). Here, we are interested in possible sources of malfunction, i.e. output values which can lead to performance
degradation or even to the failure of the IC. The maximum temperature Tmax
developed in the IC’s bondwires is the most probable cause for such malfunctions, hence we choose Tmax to be the QoI. Temperature Tmax is a functional
of the discrete solution of the electrothermal problem, i.e. Tmax = f (T). A
bondwire’s temperature depends on the random relative elongation described
above through (12) and (13). Hence, denoting a QoI with Q, in the present
setting it holds that
Q(Y) = Tmax (Y) = f (T(Y)) .

4

(14)

High Dimensional Uncertainty Quantification

We are interested in computing statistical moments of a QoI Q(Y), such as its
expectation value E[Q] = µQ and its standard deviation σQ . A common way to
6

compute these is with sampling methods, e.g. the Monte Carlo (MC) method,
which however converges slowly. While more efficient methods are available,
they suffer from the so-called curse of dimensionality, i.e. their computational
work increases, often exponentially, with respect to the number of input RVs.
We are therefore interested in methods that either alleviate or break the curse,
while maintaining good convergence properties. Moreover, we are interested
in non-intrusive methods, i.e. methods that allow the use of deterministic numerical methods for the solution of the semi-discrete problem, i.e. after its
discretization in space and time, without any modifications. Two methods have
been found to fulfill these criteria and are presented in the following sections,
namely the stochastic collocation (SC) on sparse grids (SG) and the tensor train
(TT) decomposition.

4.1

Stochastic Collocation on Sparse Grids

In the SC method, introduced in [2, 26], the deterministic problem is solved
exactly on a set of RV realizations which are called collocation points. The
choice of the collocation points y(i) is based on the distributions of the input
RVs. We seek to construct a global polynomial approximation of the QoI Q,
such that
M
X
Q(y) ≈ Q̃(y) =
αm Lm (y),
(15)
m=1

for suitable multivariate polynomials {Lm (y)}M
m=1 and coefficients αm . A brief
summary of the method follows, for more extensive presentations see [2, 17].
First, we introduce a tensor-product space P` , which is the span of tensorproduct polynomials with degree at most p = [p1 , . . . , pN ]. The positive integer
NN
` is called the approximation level. It holds that P` = n=1 Ppn , where Ppn
(i )
is the span of univariate polynomials up to degree pn . We denote with yn n ,
0 ≤ in ≤ pn , the (pn + 1) roots of the polynomials of degree (pn + 1) living in
Ppn , which are orthogonal with respect to the univariate PDF ρn (yn ).
Next, let a multi-index j = [j1 , . . . , jN ], such that the polynomial degrees pn
are given as functions of the indices jn , i.e. pn = p(jn ). The multi-dimensional
Lagrange interpolation operator corresponding to the multi-index j reads
p(j1 )

Ij,N Q (y) =

X
i1 =0

p(jN )

···

X

N
Y

(i )
(i )
Q y1 1 , . . . , yNN
ln,in (yn ) ,

(16)

n=1

iN =0

p(j )

n
where {ln,in (yn )}in =0
are univariate Lagrange polynomials of degree p(jn ), such
that
p(jn )
(k)
Y
yn − yn
.
(17)
ln,in (yn ) =
(in )
(k)
− yn
k=0,k6=in yn

For p(jn ) = jn and j = [`, `, . . . , `], we obtain the isotropic full tensorproduct space, i.e. in each direction we take all polynomials of degree at most
7

N

p(`) = `. In this case, the dimension
of P` is dim [P

 ` ] = (` + 1) . Accordingly,
(i1 ) (i2 )
(iN )
(i)
the collocation points y = y1 , y2 , . . . , yN
are given as the tensor grid
(0)

(`)

(0)

(`)

(0)

(`)

TP
Hj,N
:= {y1 , . . . , y1 } × {y2 , . . . , y2 } × · · · × {yN , . . . , yN },

(18)

and the corresponding Lagrange interpolation operator is denoted with I`,N .
The global index i is associated with the local indices in as in [17].
For the computation of statistical moments, we use a multi-dimensional
Gauss quadrature formula for integral approximation. The number of quadrature points per dimension is equal to (` + 1). Then, taking as collocation
points the tensor grid created by the one-dimensional quadrature abscissas, we
write
Z
I
X
E[Q] =
Q(y)ρ(y)dy ≈
wi Q(y(i) ),
(19)
Ω

i=1

where Ω denotes the set of random outcomes and the weights wi are given as
wi =

N
Y

Z
win

and

win =
Ωn

n=1

2
ln,i
(yn )ρn (yn )dyn .
n

However, the tensor grid approach soon becomes intractable, as dim [P` ] and
the number of collocation points increase exponentially with the number of RVs.
Therefore, collocation in high dimensions exploits sparse grids [3, 6, 18]. For
simplicity, only isotropic total degree and Smolyak grids are presented. Readers
interested in anisotropic grids are referred to [19] and the references therein.
For the total degree case, we set p(j) = j, while for the Smolyak grid case we
set p(j) = 2j for j > 0, and p(0) = 0. For both cases we require that |j| ≤ p(`),
PN
where |j| = n=1 jn . Then, the sparse interpolation is obtained by the isotropic
Smolyak formula


X
N −1
A`,N Q (y) =
(−1)`−|j|
Ij,N Q (y) .
(20)
` − |j|
`−N +1≤|j|≤`

TD
SM
, for the total degree
The associated sparse grids are denoted with H`,N
and H`,N
and the Smolyak case, respectively.

4.2

Tensor Train Decomposition

Considering again the full tensor-product space introduced in Section 4.1 and
the quadrature formula (19), we define the N -dimensional arrays, or N -tensors,
Q, W ∈ RI1 ×···×IN , with elements


(i )
(i ) (i )
Q (i1 , i2 , . . . , iN ) := Q y1 1 , y2 2 , · · · , yNN
and W (i1 , i2 , . . . , iN ) := wi1 wi2 · · · wiN ,
(21)
where in = 1, . . . , In . In the isotropic case, it holds that I1 = · · · = IN = (` + 1).
8

Then, (19) can be re-written in the form of the inner product between Q
and W, i.e.
X
E[Q] ≈ E[Q̃] = hQ, Wi :=
Q (i1 , . . . , iN ) W (i1 , . . . , iN ) .
(22)
i1 ,...,iN

For the comfort of the reader, the standard notation for tensor and matrix
elements is substituted by a MATLAB-like notation, such that Q(i1 , . . . , iN ) =
qi1 ,...,iN and W(i1 , . . . , iN ) = wi1 ,...,iN . Obviously, the storage needs of Q and
W and the computational work involving these tensors increases exponentially
with respect to N . However, this curse of dimensionality can be overcome
using low-rank tensor decomposition methods, which yield representations or
approximations of tensors, in formats requiring much less storage and allowing
efficient computations.
Low-rank decomposition is a problem which is well-studied for the 2-dimensional,
i.e. matrix, case. Using element-wise notation, a rank-R decomposition of a
matrix A ∈ RI1 ×I2 reads
A(i1 , i2 ) ≈ Ã(i1 , i2 ) =

R
X

U(i1 , r)V(r, i2 ).

(23)

r=1

If R < rank (A), the decomposition yields a rank-R approximation of the matrix. The best rank-R approximation in the Frobenius-norm sense can be obtained by the truncated singular value decomposition (SVD), such that Ã =
∗
, which however is expensive to compute. A computationally less exUR SR VR
pensive alternative to the SVD is the skeleton or cross approximation [4, 11],
defined as

−1

A ≈ Ã = A :, I 2 A I 1 , I 2
A I 1, : ,
(24)
or equivalently, in element-wise notation
A(i1 , i2 ) ≈ Ã(i1 , i2 ) =

R
X


−1

A i1 , i2t A i1s , i2t
A i1s , i2 ,

(25)

s,t=1
2
2 R
where I 1 = {i1s }R
s=1 and I = {it }t=1 denote row and column index sets, respectively. The approximation matrix Ã can be computed by touching only
(I1 + I2 )R − R2 entries of the original matrix. Therefore, this approach is particularly efficient for function-related matrices, where A(i1 , i2 ) := f (xi1 , xi2 ),
with x : ij ∈ N 7→ xij ∈ R and assuming that f (x1 , x2 ) is a smooth function. The accuracy of the approximation depends on the choice of index sets
1
2
I 1 , I 2 . For
 the best possible approximation, I and I must be chosen such that
1
2
A I , I is a submatrix of maximal volume, meaning that it has the maximal
in modulus determinant over all R × R submatrices of A. This task however
is computationally challenging and often unfeasible. In practice only quasi maximal volume submatrices [12] are required. The choice of index sets such
that the
or quasi-maximality criterion is satisfied, is often denoted
 maximality

with I 1 , I 2 = maxvol [A].

9

Various efforts have been made to generalize the low-rank decomposition
(23) from matrices to tensors [13, 14, 15, 16]. Here, we consider only the tensor
train (TT) decomposition [20], which for an N -tensor A ∈ RI1 ×···×IN is defined
as
A(i1 , . . . , iN ) ≈ ATT (i1 , . . . , iN ) =

R0 X
,...,RN

G1 (r0 , i1 , r1 )G2 (r1 , i2 , r2 ) · · · GN (rN −1 , iN , rN )

r0 ,...,rN

=

R Y
N
X
r

Gn (rn−1 , in , rn ) ,

(26)

n=1

where Gn ∈ RRn−1 ×In ×Rn are the TT-cores, R = (R0 , R1 , . . . , RN ) are the TTranks, in = 1, . . . , In are the mode indices and rn = 1, . . . , Rn are auxiliary
indices. It holds that R0 = RN = 1. The TT decomposition is also known as
the Matrix Product State (MPS) format, due to fact that it can be written as
ATT (i1 , . . . , iN ) = G1 (i1 )G2 (i2 ) · · · GN (iN ),

(27)

where Gn (in ) are Rn−1 × Rn matrices, one for each index in .
We define the multi-indices i≤n := (i1 , . . . , in ) and i>n := (in+1 , . . . , iN ),
and the n-th unfolding or matricization as the matrix An ∈ R(I1 ···In )×(In+1 ···IN )
with elements

(28)
An (i≤n , i>n ) = An i1 , . . . , in , in+1 , . . . , iN = A(i1 , . . . , iN ).
The TT decomposition is exact for Rn = rank (An ), and approximate for Rn <
rank (An ). Setting max(In ) = I and max(Rn ) = R, n = 1, . . . , N , the storage
of a tensor in the TT format requires O N IR2 elements, which is substantially
less than the O I N entries of the full tensor and especially helpful in the case
of function-related tensors. Storage and computational operations in the TT
format scale linearly with the number of dimensions N , thus breaking the curse
of dimensionality [20].
The TT approximation is typically computed by solving the minimization
problem
min kA − ATT kF ,

G1 ,...,GN

(29)

where k · kF denotes the Frobenius norm. In [21], an approach based on the
alternating least squares (ALS) algorithm is proposed to solve (29). The idea
relies on organizing the ALS iteration in left-to-right and right-to-left sweeps
and fixing at each step all TT-cores but one, thus reducing the number of
unknowns to Rn−1 × In × Rn at a time. In this, so-called, TT-cross algorithm,
cross approximation replaces least-squares in optimizing the TT-core at each
Rn
Rn
>n
step. We denote with I ≤n = {i≤n
= {i>n
tn }tn =1 the index sets
sn }sn =1 and I
with the positions of the Rn rows and columns of matricization An , such that






>n Rn
= A I ≤n , I >n . Further assuming that I ≤n , I >n =
A i≤n
sn , itn
s ,t =1
n

n

10

maxvol [An ], the generalized TT-cross approximation from [21] reads
ATT (i1 , . . . , iN ) =

N
XY




>n
>n −1
A is≤k−1
,
i
,
i
A i≤n
.
n
tn
sn , itn
n−1

(30)

s,t n=1

Border sets I ≤0 = ∅ and I >N = ∅ have been introduced tounify the notation.
The total complexity of the TT-cross algorithm is O N IR3 , i.e. scales linearly
with respect to the number of dimensions N .
TT-cross has the drawback that the TT-ranks Rn have to be guessed in
advance. Underestimations of the TT-ranks may lead to poor approximations,
while overestimations increase computational costs. A more efficient and rankrevealing method is presented in [22]. The corresponding algorithm is based on
the Density Matrix Renormalization Group (DMRG) numerical technique from
quantum chemistry, which can be seen as an idea of constrained optimization
in the MPS format. In this TT-RC (Renormalization Cross) algorithm, optimization is performed over two TT-cores, Gn and Gn+1 , simultaneously. At step
n of the sweeps, the supercore Wn (in , in+1 ) = Gn (in ) Gn+1 (in+1 ) is found by
solving (29). Matrix Wn is of size Rn−1 In ×In+1 Rn+1 , hence it does not include
any information regarding rank Rn . Matrices Gn and Gn+1 are recovered using
a truncated-SVD of Wn . This, so-called, decimation step provides us with the
TT-rank Rn , thus the rank-revealing property of the method.
In [23], a greedy approach to TT-cross approximation is discussed, similar
to the adaptive cross approximation (ACA) [4] method for matrices. The proposed algorithm, called greedy-TT-cross, modifies the TT-RC algorithm and
adds at each step a new pivot index, based on where the error of the current
approximation is quasi-maximum in the current step’s supercore Wn . Its total
computational complexity is also O(N IR3 ), i.e. it remains linear with respect
to the number of dimensions N .

5

Numerical Results

Simulations of the electrothermal field problem presented in Section 2 were
performed. The implicit Euler method was used for time integration. In total,
the initial 3.5 seconds of the IC’s duty cycle were simulated, divided into 51
equidistant time-steps. The total amount of degrees of freedom due to the FIT
discretization was 10516, for each time-step.
First, a Monte Carlo (MC) simulation with 40000 realizations of the input
RVs was performed. Each realization corresponds to a single call of the deterministic solver which solves (7) and (8). For both E [Tmax ] and σTmax we have
convergence up to the second digit. The most accurate values of the mean and
the standard deviation of Tmax , i.e. those computed with the full set of 40000
MC
MC
samples, are E [Tmax ]
= µMC
Tmax = 544.42 K and σTmax = 4.07 K, respectively.
The main question is whether we can reach the most accurate MC mean and
standard deviation values with significantly less calls to the solver, by applying
either the SC-SG or the TT decomposition.
11

The SC-SG method was implemented using the Dakota toolkit [27]. Approximation levels ` = 1, 2, 3 have been considered for the case of Smolyak sparse
grids. The quadrature abscissas chosen for (19) are the Clenshaw-Curtis (CC).
The number of calls to the deterministic solver is equal to the number of collocation points. The relative errors for the mean value and the standard deviation
MC
are computed with µMC
Tmax and σTmax as the reference values, respectively. The
corresponding data and results for all cases are presented in Table 1.
`
1
2
3

Collocation points
25
337
3249

St. dev. rel. error σ,SG (%)
3.50
3.19
3.36

Table 1: Results for the stochastic collocation on sparse grids.
The relative errors for the mean value are consistently below 1% and therefore omitted. The relative errors for the standard deviation are consistently
between 3% − 4%, even for a the case of ` = 1, i.e. with only 25 collocation
points. Therefore, a good approximation of the reference mean and standard
deviation values can be obtained using stochastic collocation on sparse grids,
with significantly less solver calls. However, some oscillation in the relative error
of the standard deviation still remains, therefore it can be deduced that more
approximation levels are required until full convergence is reached.
For the implementation of the TT decomposition approach, the MATLAB
TT-Toolbox [28] was used. We denote with ⊗ the outer product and with
wn ∈ R(`+1) the vectors containing all quadrature weights for the corresponding
RVs Yn , such that wn (in ) = win . Then, the N -tensor W from (21) can be
written as
W = w1 ⊗ · · · ⊗ wN ,
(31)
or equivalently, in element-wise notation
W (i1 , . . . , iN ) = w1 (i1 ) · · · wN (iN ).

(32)

The formats in (31), (32) are equivalent to an exact TT decomposition, i.e.
WTT = W, with all TT-ranks equal to 1. Therefore, storing WTT requires
only N (` + 1) elements.
Since WTT is an exact TT decomposition, the computation of statistical
moments relies predominantly on the approximation of the N -tensor Q, which
is obtained using the greedy-TT-cross algorithm. The algorithm is run for an
increasing number of allowed sweeps, which in turn leads to increasing TT-ranks,
larger interpolation index sets and thus better approximations, at the cost of
more calls to the deterministic solver. Each approximation QTT yields different
results for E[Tmax ] = µTmax and σTmax . The respective results are presented in
Table 2.
As can be observed, the maximum TT-rank increases with the number of
allowed sweeps. The maximum TT-rank is an indicator of the storage needs for
12

Sweeps
–
1
2
3
4
5
6
7
8
9
10

Solver calls
`=1 `=2 `=3
154
242
327
291
491
670
460
772
1133
647
1130 1616
846
1529 2200
1018 1938 2980
1264 2790 3713
1497 3284 4536
1753 3981 5833
2028 4550 6875

Max. TT-rank Rn,max
–
2
3
4
5
6
7
8
9
10
11

St. dev. rel. error σ,TT (%)
`=1
`=2
`=3
1.96
12.25
18.86
6.52
2.56
14.03
6.90
125.43
13.63
56.42
11.62
10.66
3.05
1.49
12.49
1.05
1.49
5.30
< 1.0
2.70
4.25
< 1.0 < 1.0
2.26
1.10
2.20
3.19
< 1.0
1.48
1.90

Table 2: Results for the tensor train decomposition.

2
Q in the TT format. As discussed in Section 4.2, O N (` + 1) Rn,max
elements
are required to store a tensor in the TT format. Here, storing QTT requires
much less space than this limit, because the maximum TT-rank appears only
in few cores.
As in the sparse grids case, the mean value relative errors are consistently
below 1% for all considered approximation levels and allowed sweeps, and are
therefore omitted. Observing the standard deviation relative errors, it seems
that accurate approximations of tensor Q are obtained when the allowed sweeps
are more than 5 for ` = 1, 2 and more than 7 for ` = 3. After that point,
deviations from the reference value are consistently less than 5%. Hence, after
a certain number of sweeps, the algorithm has located interpolation index sets
which provide a good TT approximation QTT ≈ Q, thus leading to accurate
mean and standard deviation values. It can be observed that accuracies superior
to those of the SC-SG can be obtained, albeit at the cost of more solver calls.
It can be observed that, in the TT-based approach, the most accurate results are obtained for ` = 1, which indicates low-level interactions between the
RVs. This should be investigated further, since it could be the case that high
approximation levels do not really contribute to better accuracies. A possible
question is whether the use of full tensor-product spaces of lower order would
yield better approximation results. With respect to the more accurate results of
the TT decomposition approach, it should be noted that isotropic sparse grids
are a priori constructed for a class of smooth-enough functions, while the TT
approximation is built for a specific function. This is a structural disadvantage
for the SC-SG method, which should be taken into consideration. Adaptive construction of sparse grids, i.e. grids tailored to a specific function, is a possible
extension of this work.

13

6

Conclusion

Uncertainty quantification methods based on stochastic collocation on sparse
grids and on tensor train decompositions were implemented for an electrothermal problem with 12 random variables. Accurate results for the first two
statistical moments, i.e. mean and standard deviation, of the considered quantity of interest were obtained with both methods, with significantly less solver
calls than using the Monte Carlo method. Stochastic collocation on sparse grids
appears to be the computationally most efficient method for the considered
problem setting. At the cost of higher computational burden, the tensor train
approach yields results of superior accuracy, compared to stochastic collocation.
However, the aforementioned structural advantage of the tensor train decomposition has to be taken into consideration. Extensions of the present work
shall handle problems involving an increased number of random variables, with
the goal of determining the efficiency of both approaches, especially of the more
novel tensor train decomposition, when the impact of the curse of dimensionality
is more pronounced.

Acknowledgement
This work was supported by the Graduate School of Computational Engineering,
Technische Universität Darmstadt. Part of this work was supported by the fp7nanoCOPS (Nanoelectronic Coupled Problems Solutions) project.

References
[1] Babuška I., Tempone R., and Zouraris G.E. 2004. Galerkin finite element
approximations of stochastic elliptic partial differential equations. SIAM J.
Numer. Anal. 42, 800-825.
[2] Babuška I., Nobile F. and Tempone R. 2010. A stochastic collocation method
for elliptic partial differential equations with random input data. SIAM Rev.
52, 317-355.
[3] Barthelmann V., Novak E. and Ritter K. 2000. High dimensional polynomial interpolation on sparse grids. Advances in Computational Mathematics 12, 273-288.
[4] Bebendorf M. 2000. Approximation of boundary element matrices. Numer.
Math. 86, 565-589.
[5] Bellman R. 1957. Dynamic programming princeton university press. Princeton, NJ.
[6] Bungartz H.J. and Griebel M. 2004. Sparse grids. Acta numerica 13, 147269.

14

[7] Nguyen T.T., Mac D.H., and Clénet S. 2015. Non intrusive uncertainty quantification method for models with a high number of parametersApplication to a magnetoelectric sensor. Proceedings of Conference on the
Computation of Electromagnetic Fields (COMPUMAG).
[8] Calisch R.E. 1998. Monte Carlo and Quasi Monte Carlo methods. Acta
Numerica 7, 1-49.
[9] Casper T., Schöps S. and De Gersem H. 2016. Electrothermal simulation
of bonding wire degradation under uncertain geometries. Proceedings of
Conference on Design, Automation and Test in Europe (DATE).
[10] Chkifa A., Cohen A., and Schwab C. 2015. Breaking the curse of dimensionality in sparse polynomial approximation of parametric PDEs. Journal
de Mathématiques Pures et Appliquées 103, 400-428.
[11] Goreinov S.A., Tyrtyrshnikov E.E., Zamarashkin N.L. 1997. A theory of
pseudo-skeleton approximations. Linear Algebra Appl. 261, 1-21.
[12] Goreinov S.A., Oseledets I.V., Savostyanov D.V., Tyrtyrshnikov E.E. and
Zamarashkin N.L. 2010. How to find a good submatrix. Matrix Methods:
Theory, Algorithms, Applications, World Scientific, Hackensack NY, 247256.
[13] Grasedyck L., Kressner D. and Tobler C. 2013. A literature survey of lowrank tensor approximation techniques, GAMM-Mitteilungen 36(1), 53-78.
[14] Hackbusch W. 2012. Tensor spaces and numerical tensor calculus. SpringerVerlag, Berlin.
[15] Khoromskij B.N. 2012. Tensor-structured numerical methods in scientific
computing: survey on recent advances. Chemometr. Intell. Lab. Syst. 110,
1-19.
[16] Kolda T.G. and Bader B.W. 2009. Tensor decompositions and applications.
SIAM Rev. 51, 455-500.
[17] Motamed M., Nobile F. and Tempone R. 2013. A stochastic collocation
method for the second order wave equation with a discontinuous random
speed. Numer. Math. 123, 493-536.
[18] Nobile F., Tempone R. and Webster C.G. 2008. A sparse grid stochastic
collocation method for partial differential equations with random input data.
SIAM J. Numer. Anal. 46, 23092345.
[19] Nobile F., Tempone R. and Webster C.G. 2008. An anisotropic sparse grid
stochastic collocation method for partial differential equations with random
input data. SIAM J. Numer. Anal. 46, 2411-2442.
[20] Oseledets I.V. 2011. Tensor-train decomposition. SIAM J. Sci. Comput. 33,
2295-2317.
15

[21] Oseledets I.V. and Tyrtyrshnikov E.E. 2010. TT-cross approximation for
multidimensional arrays. Linear Algebra Appl. 432, 70-88.
[22] Savostyanov D.V., Oseledets I.V. 2011. Fast adaptive interpolation of multidimensional arrays in tensor train format. Proceedings of IEEE 7th International Workshop on Multidimensional Systems (nDS).
[23] Savostyanov D.V. 2014. Quasioptimality of maximum-volume cross interpolation of tensors. Linear Algebra Appl. 458, 217-244.
[24] Uschmajew A. 2013. Zur Theorie der Niedrigrangapproximation in Tensorprodukten von Hilberträumen. PhD Thesis, Technische Universität Berlin.
[25] Weiland T. 1996. Time domain electromagnetic field computation with finite
difference methods. Int. J. Numer. Model. El. 9, 295-319.
[26] Xiu D. and Hesthaven J.S. 2005. High-order collocation methods for differential equations with random inputs. SIAM J. on Sci. Comput. 27, 11181139.
[27] Adams B.M., Bauman L.E., Bohnhoff W.J., Dalbey K.R., Ebeida M.S.,
Eddy J.P., Eldred M.S., Hough P.D., Hu K.T., Jakeman J.D., Stephens J.A., Swiler L.P., Vigil D.M., and Wildey T.M. 2014. Dakota, A Multilevel Parallel Object-Oriented Framework for Design Optimization, Parameter Estimation, Uncertainty Quantification, and Sensitivity Analysis:
Version 6.0 User’s Manual. Sandia Technical Report SAND2014-4633. Updated November 2015 (Version 6.3).
[28] TT-toolbox. https://github.com/oseledets/TT-Toolbox.

16

