[PPT] - Uncertainty Propagation in Linear Systems: An Exact Solution Using PowerPoint Presentation

SLIDE 1

Waikiki, Hawaii, 24 April 2007

Uncertainty Propagation in Linear Systems: An Exact Solution Using random Matrix Theory

S Adhikari

School of Engineering, University of Wales Swansea, Swansea, U.K. Email: S.Adhikari@swansea.ac.uk URL: http://engweb.swan.ac.uk/∼adhikaris

Inverse of a Random Matrix – p.1/36

SLIDE 2

Waikiki, Hawaii, 24 April 2007

Outline

Motivation Current methods for response-statistics calculation Matrix variate probability density functions Exact inverse of a general real symmetric random matrix Exact response moments of linear systems Numerical example Conclusions

Inverse of a Random Matrix – p.2/36

SLIDE 3

Waikiki, Hawaii, 24 April 2007

Background

In many stochastic mechanics problems we need to solve a system of linear stochastic equations: Ku = f. (1) K ∈ Rn×n is a n × n real non-negative definite random matrix, f ∈ Rn is a n-dimensional real deterministic input vector and u ∈ Rn is a n-dimensional real uncertain

utput vector which we want to determine.

This typically arise due to the discretisation of stochastic partial differential equations (eg. in the stochastic finite element method)

Inverse of a Random Matrix – p.3/36

SLIDE 4

Waikiki, Hawaii, 24 April 2007

Background

In the context of linear structural mechanics, K is known as the stiffness matrix, f is the forcing vector and u is the vector of structural displacements. Often, the objective is to determine the probability density function (pdf) and consequently the cumulative distribution function (cdf) of u. This will allow one to calculate the reliability of the system. It is generally difficult to obtain the probably density function (pdf) of the response. As a consequence, engineers often intend to obtain only the fist few moments (typically the fist two) of the response quantity.

Inverse of a Random Matrix – p.4/36

SLIDE 5

Waikiki, Hawaii, 24 April 2007

Objectives

We propose an exact analytical method for the inverse of a real symmetric (in general non-Gaussian) random matrix of arbitrary dimension. The method is based on random matrix theory and utilizes the Jacobian of the underlying nonlinear matrix transformation. Exact expressions for the mean and covariance of the response vector is obtained in closed-form.

Inverse of a Random Matrix – p.5/36

SLIDE 6

Waikiki, Hawaii, 24 April 2007

Current Approaches

The random matrix can be represented as K = K0 + ∆K (2) K0 ∈ Rn×n is the deterministic part and the random part: ∆K =

m

j=1

ξjKI

j + m

j=1

m

l=1

ξjξlKII

jl + · · ·

(3) m is the number of random variables, KI

j, KII jl ∈ Rn×n, ∀ j, l

are deterministic matrices and ξj, ∀ j are real random variables.

Inverse of a Random Matrix – p.6/36

SLIDE 7

Waikiki, Hawaii, 24 April 2007

Perturbation based approach

Represent the response as u = u0 + ξjuI

j + m

j=1

m

l=1

ξjξluII

jl + · · · .

(4) where u0 = K0−1f (5) uI

j = −K0−1KI ju0,

∀ j (6) and uII

jl = −K0−1[KII jl u0 + KI juI l + KI l uI j],

∀ j, l. (7)

Inverse of a Random Matrix – p.7/36

SLIDE 8

Waikiki, Hawaii, 24 April 2007

Neumann expansion

Provided

K0−1∆K
F < 1,

K−1 =

K0(In + K0−1∆K)

−1 = K0−1 − K0−1∆KK0−1 + K0−1∆KK0−1∆KK0−1 + · · · . Therefore, u = K−1f = u0 − Tu0 + T2u0 + · · · (8) where T = K0−1∆K ∈ Rn×n is a random matrix.

Inverse of a Random Matrix – p.8/36

SLIDE 9

Waikiki, Hawaii, 24 April 2007

Projection methods

Here one ‘projects’ the solution vector onto a complete stochastic basis. Depending on how the basis is selected, several methods are proposed. Using the classical Polynomial Chaos (PC) projection scheme u =

P−1

j=0

ujΨj(ξ) (9) where uj ∈ Rn, ∀j are unknown vectors and Ψj(ξ) are multidimensional Hermite polynomials in ξr.

Inverse of a Random Matrix – p.9/36

SLIDE 10

Waikiki, Hawaii, 24 April 2007

A partial summary

Methods Sub-methods 1. Perturbation First and second order perturbation 1,2, based methods Neumann expansion 3,4, improved perturbation method 5. 2. Projection methods Polynomial chaos expansion 6, random eigenfunction expansion 4, stochastic reduced basis method 7–9, Wiener−Askey chaos expansion 10–12, domain decomposition method 13,14. 3. Monte carlo simulation Simulation methods 15,16, and other methods Analytical method in references 17–21, Exact solutions for beams 22,23.

Inverse of a Random Matrix – p.10/36

SLIDE 11

Waikiki, Hawaii, 24 April 2007

Matrix variate distributions

The probability density function of a random matrix can be defined in a manner similar to that of a random variable. If A is an n × m real random matrix, the matrix variate probability density function of A ∈ Rn,m, denoted as pA(A), is a mapping from the space of n × m real matrices to the real line, i.e., pA(A) : Rn,m → R.

Inverse of a Random Matrix – p.11/36

SLIDE 12

Waikiki, Hawaii, 24 April 2007

Gaussian random matrix

The random matrix X ∈ Rn,p is said to have a matrix variate Gaussian distribution with mean matrix M ∈ Rn,p and covariance matrix Σ ⊗ Ψ, where Σ ∈ R+

n and Ψ ∈ R+ p provided

the pdf of X is given by pX (X) = (2π)−np/2 |Σ|−p/2 |Ψ|−n/2 etr

−1

2Σ−1(X − M)Ψ−1(X − M)T

(10)

This distribution is usually denoted as X ∼ Nn,p (M, Σ ⊗ Ψ).

Inverse of a Random Matrix – p.12/36

SLIDE 13

Waikiki, Hawaii, 24 April 2007

Symmetric Gaussian matrix

If Y ∈ Rn×n is a symmetric Gaussian random matrix then its pdf is given by pY (Y) = (2π)−n(n+1)/4 BT

n(Σ ⊗ Ψ)Bn

−1/2

etr

−1

2Σ−1(Y − M)Ψ−1(Y − M)T

.

(11) This is denoted as Y = YT ∼ SNn,n

M, BT

n(Σ ⊗ Ψ)Bn

. The

elements of the translation matrix Bn ∈ Rn2×n(n+1)/2 are: (Bn)ij,gh = 1 2 (δigδjh + δihδjg) , i ≤ n, j ≤ n, g ≤ h ≤ n, (12)

Inverse of a Random Matrix – p.13/36

SLIDE 14

Waikiki, Hawaii, 24 April 2007

Matrix variate Gamma distribution

A n × n symmetric positive definite matrix random W is said to have a matrix variate gamma distribution with parameters a and Ψ ∈ R+

n , if its pdf is given by

pW (W) =

Γn (a) |Ψ|−a−1 |W|a− 1

2(n+1) etr {−ΨW} ; ℜ(a) > 1

2(n−1) This distribution is usually denoted as W ∼ Gn(a, Ψ). Here the multivariate gamma function: Γn (a) = π

1 4n(n−1)

n

k=1

Γ

a − 1

2(k − 1)

; for ℜ(a) > (n − 1)/2

Inverse of a Random Matrix – p.14/36

SLIDE 15

Waikiki, Hawaii, 24 April 2007

Wishart matrix

A n × n symmetric positive definite random matrix S is said to have a Wishart distribution with parameters p ≥ n and Σ ∈ R+

n , if its pdf is given by

pS (S) =

2

1 2np Γn

1 2p

|Σ|

1 2p

−1 |S|

1 2(p−n−1)etr

−1

2Σ−1S

(13)

This distribution is usually denoted as S ∼ Wn(p, Σ). Note: Gn(a, Ψ) = Wn(2a, Ψ−1/2), so that Gamma and Wishart are equivalent distributions.

Inverse of a Random Matrix – p.15/36

SLIDE 16

Waikiki, Hawaii, 24 April 2007

Inverse of a scalar

ku = f (14) where k, u, f ∈ R. Suppose the pdf k is pk(k) and we are interested in deriving the pdf of h = k−1. (15) The Jacobian of the above transformation J = ˛ ˛ ˛ ˛ ∂h ∂k ˛ ˛ ˛ ˛ = ˛ ˛−k−2˛ ˛ = |k|−2 . (16) Using the Jacobian, the pdf of h can be obtained as ph(h)(dh) = pk(k)(dk) (17)

r

ph(h) = 1 ˛ ˛ ˛ ∂h

∂k

˛ ˛ ˛ pk(k) (18)

r

ph(h) = 1 J (k = h−1)pk ` k = h−1´ = |h|−2 pk ` h−1´ . (19)

Inverse of a Random Matrix – p.16/36

SLIDE 17

Waikiki, Hawaii, 24 April 2007

The case of n × n matrices

Suppose the matrix variate probability density function of the non-singular matrix K is given by pK (K) : Rn×n → R. Our interest is in the pdf (i.e joint pdf of the elements) of H = K−1 ∈ Rn×n. (20) The elements of H are complicated non-linear function of the elements of K (i.e. even if the elements of K are joint Gaussian, the elements of H will not be joint Gaussian). H may not have any banded structure even if K is of banded nature.

Inverse of a Random Matrix – p.17/36

SLIDE 18

Waikiki, Hawaii, 24 April 2007

Pdf transforation in matrix space

The procedure to obtain the pdf of H is very similar to that of the univariate case: pH (H) (dH) = pK(K)(dK) (21)

r

pH (H) = 1

dH

dK

pK(K)

(22)

r

pH (H) = 1 J

K = H−1pK
K = H−1

(23) = |H|−(n+1) pK

H−1

. (24) For the univariate case (n = 1) Eq. (24) reduces to the familiar equivalent expression obtained in Eq. (19).

Inverse of a Random Matrix – p.18/36

SLIDE 19

Waikiki, Hawaii, 24 April 2007

Derivation of the Jacobian - 1

We have KK−1 = KH = In. (25) Taking the matrix differential (dK) H + K (dH) = On

r

(dH) = −K−1 (dK) K−1. (26) Treat (dH) , (dK) ∈ Rn×n as variables and K as constant since it does not contain (dH) or (dK). Taking the vec of Eq. (26) vec (dH) = −vec

K−1 (dK) K−1

= −

K−1 ⊗ K−1

vec (dK) . (27)

Inverse of a Random Matrix – p.19/36

SLIDE 20

Waikiki, Hawaii, 24 April 2007

Derivation of the Jacobian - 2

Because (dH) and (dK) are symmetric matrices we need to eliminate the ‘duplicate’ variables appearing in the preceding linear transformation. This can be achieved in a systematic manner by using the translation matrix Bn as vecp (dH) = B†

nvec (dH) = −

B†

n

K−1 ⊗ K−1

Bn

vecp (dK) .

(28) The Jacobian associated with the above linear transformation is simply the determinant of the matrix B†

n

K−1 ⊗ K−1

Bn, that is J =

B†

n

K−1 ⊗ K−1

Bn

= |K|−(n+1) .

(29)

Inverse of a Random Matrix – p.20/36

SLIDE 21

Waikiki, Hawaii, 24 April 2007

RM model for stiffness matrix

If the mean of K is K, then K ∼ Wn (p, Σ), where p = n + 1 + θ Σ = K/α θ = 1 δ2

K

1 + {Trace
K
}2/Trace
K

2

− (n + 1) and α =

θ(n + 1 + θ).

δK is the normalized standard-deviation of K: δ2

K = E

K − E [K] 2

F

E [K] 2

F

. (30)

Inverse of a Random Matrix – p.21/36

SLIDE 22

Waikiki, Hawaii, 24 April 2007

Pdf of K−1

The pdf H = K−1, that is, the joint pdf of all the elements of H can be obtained as pH (H) = |H|−(n+1) pK

H−1

= |H|−(n+1)

2

1 2np Γn

1 2p

|Σ|

1 2 p

−1 H−1

1

2(p−n−1) etr

−1

2Σ−1H−1

= |H|−(n+1+p)/2
2

1 2np Γn

1 2p

|Σ|

1 2p

−1 etr

−1

2Σ−1H−1

.

(31) Using this exact pdf, the moments of the inverse matrix can be

btained.

Inverse of a Random Matrix – p.22/36

SLIDE 23

Waikiki, Hawaii, 24 April 2007

Inverted Wishart matrix

A n × n symmetric positive definite random matrix V is said to have an inverted Wishart distribution with parameters m and Ψ ∈ R+

n , if its pdf is given by

pV (V) = 2− 1

2 (m−n−1)n|Ψ| 1 2 (m−n−1)

Γn 1

2(m − n − 1)

|V|m/2etr
−1

2V−1Ψ

; m > 2n, Ψ > 0.

(32) This distribution is usually denoted as V ∼ IWn(m, Ψ). We can show that K−1 ∼ IWn(θ + 2n + 2, αK

−1).

Inverse of a Random Matrix – p.23/36

SLIDE 24

Waikiki, Hawaii, 24 April 2007

Moments of K−1

The first moment (mean), second-moment and the elements

f the covariance tensor of K−1 can be obtained24 exactly in

closed-form as E

K−1

= Ψ m − 2n − 2 = α θ K

−1

(33) E

K−2

= Trace (Ψ) Ψ + (m − 2n − 2)Ψ2 (m − 2n − 1)(m − 2n − 2)(m − 2n − 4) E

K−1AK−1

= Trace (AΨ) Ψ + (m − 2n − 2)ΨAΨ (m − 2n − 1)(m − 2n − 2)(m − 2n − 4) = α2 Trace

AK

−1

K

−1 + θK −1AK −1

θ(θ + 1)(θ − 2) . (34)

Inverse of a Random Matrix – p.24/36

SLIDE 25

Waikiki, Hawaii, 24 April 2007

Response Moments - 1

The complete response vector is u = K−1f. (35) In many practical problems only few elements of u or linear combinations of some elements of u may be of

interest. Therefore, we are interested in the quantity

y = Ru = RK−1f; R ∈ Rr×n (36) The matrix R can be also selected to ‘extract’ other physical quantities such as the stress components within

ne element or a group of elements.

Inverse of a Random Matrix – p.25/36

SLIDE 26

Waikiki, Hawaii, 24 April 2007

Response Moments - 2

The mean of y: ¯ y = E [y] = E

RK−1f
= R E
K−1

f = α θ RK

−1f.

(37) The complete covariance matrix of y: cov (y, y) = E

(y − ¯

y)(y − ¯ y)T = E

yyT

− ¯ y¯ yT = R E

K−1ffTK−1

RT − ¯ y¯ yT = α2Trace

ffTK

−1

RK

−1RT + θ(θ + 2)¯

y¯ yT θ(θ + 1)(θ − 2) . (38)

Inverse of a Random Matrix – p.26/36

SLIDE 27

Waikiki, Hawaii, 24 April 2007

Steps for complete analysis

Obtain normalized standard deviation δ2

G = E h

K−E[K]

2 F

i

E[K]

2 F

=

Trace(cov(vec(K))) Trace „

K

2«

. Calculate the constants θ = 1 δ2

K

1 + {Trace
K
}2/Trace
K

2

− (n + 1), p = n + 1 + θ, α =

θ(n + 1 + θ) and Σ = K/α.

The mean: ¯ y = α θ RK

−1f.

The covariance: cov (y, y) = α2Trace

ffTK

−1

RK

−1RT + θ(θ + 2)¯

y¯ yT θ(θ + 1)(θ − 2) .

Inverse of a Random Matrix – p.27/36

SLIDE 28

Waikiki, Hawaii, 24 April 2007

Example: A cantilever Plate

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 −0.5 0.5 1

4 3 5

X direction (length)

6

Outputs

1

Input, f=10 N

2

Y direction (width) Fixed edge

A steel cantilever plate with a slot; ¯ E = 200 × 109N/m2, ¯ µ = 0.3, ¯ t = 7.5mm, Lx = 1.2m, Ly = 0.8m; 25 × 15 elements resulting n = 1200.

Inverse of a Random Matrix – p.28/36

SLIDE 29

Waikiki, Hawaii, 24 April 2007

Stochastic Properties

The Young’s modulus, Poissons ratio and thickness are random fields of the form E(x) = ¯ E (1 + ǫEf1(x)) (39) µ(x) = ¯ µ (1 + ǫµf2(x)) (40) and t(x) = ¯ t (1 + ǫtf4(x)) (41) The strength parameters are: ǫE = 0.15, ǫµ = 0.10, and ǫt = 0.15. The random fields fi(x), i = 1, 2, 3 are delta-correlated homogenous Gaussian random fields.

Inverse of a Random Matrix – p.29/36

SLIDE 30

Waikiki, Hawaii, 24 April 2007

Response calculation

The value of δk (calculated using a 5000-sample Monte Carlo simulation of the random fields) is obtained as δK = 0.2616. From the 1200 × 1200 stiffness matrix we obtain Trace

K
= 5.5225×109

and Trace

K

2

= 9.6599×1016. and θ = 3.4274 × 103 and α = 3.9827 × 103.

Inverse of a Random Matrix – p.30/36

SLIDE 31

Waikiki, Hawaii, 24 April 2007

Comparison of results

The mean and standard deviation of the response vector. The numbers in the parenthesis correspond to the percentage error in the Monte Carlo Simulation (with 1000 samples) results with respect to the exact analytical results.

Response Analytical MSC Analytical standard MSC standard quantity mean (mm) mean (mm) deviation (mm) deviation (mm) y1 = u112 5.5058 5.5178 (0.218 %) 0.1438 0.1459 (1.436 %) y2 = u325 2.6420 2.6475 (0.208 %) 0.0734 0.0740 (0.818 %) y3 = u658 10.2265 10.2485 (0.216 %) 0.2537 0.2561 (0.972 %) y4 = u1045 12.6039 12.6317 (0.221 %) 0.3294 0.3313 (0.570 %) y5 = u868 5.9608 5.9725 (0.197 %) 0.1586 0.1604 (1.155 %) y6 = u205 10.2951 10.3169 (0.212 %) 0.2507 0.2547 (1.580 %)

Inverse of a Random Matrix – p.31/36

SLIDE 32

Waikiki, Hawaii, 24 April 2007

Convergence of results - 1

0.5 1 1.5 2 2.5 3 3.5 4 x 104 0.02 0.04 0.06 0.08 0.1 0.12 0.14 Number of samples Percentage error in mean y1 y2 y3 y4 y5 y6

Variation of the relative error in MCS mean with respect to the number of samples.

Inverse of a Random Matrix – p.32/36

SLIDE 33

Waikiki, Hawaii, 24 April 2007

Convergence of results - 2

0.5 1 1.5 2 2.5 3 3.5 4 x 104 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Number of samples Percentage error in standard deviation y1 y2 y3 y4 y5 y6

Variation of the relative error in MCS standard deviation with respect to the number of samples.

Inverse of a Random Matrix – p.33/36

SLIDE 34

Waikiki, Hawaii, 24 April 2007

Summary - 1

The probabilistic characterization of the response of linear stochastic systems requires inverse of a real symmetric random matrix (an outstanding problem for more than four decades). An exact and simple closed-form expression of the joint probability density function of the elements of the inverse

f a symmetric random matrix is derived.

A matrix itself is treated like a variable, as opposed to view it as a collection of many variables. This outlook significantly simples the calculation of the Jacobian involved in the non-linear matrix transformation.

Inverse of a Random Matrix – p.34/36

SLIDE 35

Waikiki, Hawaii, 24 April 2007

Summary - 2

The random matrices considered are in general non-Gaussian and of arbitrary dimensions. Moments of the response do not require a series/perturbation/PC expansion. The numerical implementation is straight-forward and non-intrusive.

Inverse of a Random Matrix – p.35/36

SLIDE 36

Waikiki, Hawaii, 24 April 2007

Open Issues

Any real matrix pdf can be used for pK(K) and the pdf of H = K−1 can be obtained. However, obtaining the response pdf (requires further transformation) or response moments from pH(H) is a not trivial task. Selecting a matrix variate pdf to matrix data is a challenging task itself (topic of my Thursdays paper). The inverse of a complex symmetric random matrix cannot be obtained easily from the proposed formulation. As a result, it is applicable to static or undamped systems

nly and therefore is of somewhat limited applicability.

Inverse of a Random Matrix – p.36/36

SLIDE 37

References

[1] Shinozuka M., Yamazaki F. S.T. Ariaratnam, G.I. Schueller, I. El- ishakoff, eds., Stochastic structural dynamics: Progress in theory and applica-

tions. Elsevier Applied Science, London, 1998.

[2] Kleiber M., Hien T.D. The Stochastic Finite Element Method. John Wiley, Chichester, 1992. [3] Yamazaki F., Shinozuka M., Dasgupta G. ASCE Journal of Engineering Mechanics, 114(8), 1335–1354, 1988. [4] Adhikari S., Manohar C.S. International Journal for Numerical Methods in Engineering, 44(8), 1157–1178, 1999. [5] Elishakoff I., Ren Y.J., Shinozuka M. Chaos Solitons & Fractals, 5(5), 833–846, 1995. [6] Ghanem R., Spanos P. Stochastic Finite Elements: A Spectral Approach. Springer-Verlag, New York, USA, 1991. [7] Nair P.B., Keane A.J. AIAA Journal, 40(8), 1653–1664, 2002. [8] Sachdeva S.K., Nair P.B., Keane A.J. Computer Methods in Applied Me- chanics and Engineering, 195(19-22), 2371–2392, 2006. [9] —. Probabilistic Engineering Mechanics, 21(2), 182–192, 2006. [10] Xiu D.B., Karniadakis G.E. Siam Journal on Scientific Computing, 24(2), 619–644, 2002. [11] —. Journal of Computational Physics, 187(1), 137–167, 2003. [12] Wan X.L., Karniadakis G.E. Journal of Scientific Computing, 27((-3), 455– 464, 2006. [13] Sarkar A., Benabbou N., Ghanem R. 47th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics & Materials Conference. Newport, Rhode Is- land, USA, 2006. [14] —. 7th World Congress on Computational Mechanics. Los Angeles, CA, USA, 2006. [15] Hurtado J.E., Barbat A.H. Archives of Computational Methods in Engineer- ing, 5(1), 3–29, 1998. [16] Papadrakakis M., Papadopoulos V. Computer Methods in Applied Mechanics and Engineering, 134(3-4), 325–340, 1996. [17] Muscolino G., Ricciardi G., Impollonia N. Probabilistic Engineering Me- chanics, 15(2), 199–212, 2000. [18] Impollonia N., Muscolino G. Meccanica, 37(1-2), 179–192, 2002. [19] Falsone G., Impollonia N. Computer Methods in Applied Mechanics and En- gineering, 191(44), 5067–5085, 2002. [20] —. Computer Methods in Applied Mechanics and Engineering, 192(16-18), 2187–2188, 2003. [21] Impollonia N., Ricciardi G. Probabilistic Engineering Mechanics, 21(2), 171–181, 2006. [22] Elishakoff I., Ren Y.J., Shinozuka M. International Journal of Solids and Structures, 32(16), 2315–2327, 1995. [23] Elishakoff I., Ren Y.J. Large Variation Finite Element Method for Stochastic

Problems. Oxford University Press, Oxford, U.K., 2003.

[24] Gupta A., Nagar D. Matrix Variate Distributions. Monographs & Surveys in Pure & Applied Mathematics. Chapman & Hall/CRC, London, 2000.

36-1