Independent Component Analysis to enhance performances of Karhunen–Loeve expansions for non-Gaussian stochastic processes: Application to uncertain systems

Journal of Sound and Vibration 333 (2014) 5600–5613

Contents lists available at ScienceDirect

Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi

Independent Component Analysis to enhance performances of Karhunen–Loeve expansions for non-Gaussian stochastic processes: Application to uncertain systems Mohammad Khalil 1, Abhijit Sarkar n Carleton University, Department of Civil and Environmental Engineering, Mackenzie Building, Colonel By Drive, Ottawa, Ontario, Canada K1S 5B6

a r t i c l e in f o

abstract

Article history: Received 30 April 2013 Received in revised form 18 February 2014 Accepted 11 April 2014 Handling Editor: K. Worden Available online 27 June 2014

In the context of engineering systems, an essential step in uncertainty quantification is the development of accurate and efficient representation of the random input parameters. For such input parameters modeled as stochastic processes, Karhunen–Loeve expansion is a classical approach providing efficient representations using a set of uncorrelated, but generally statistically dependent random variables. The dependence structure among these random variables may be difficult to estimate statistically and is thus ignored in many practical applications. This simplifying assumption of independence may lead to considerable errors in estimating the variability in the system state, thus limiting the effectiveness of Karhunen–Loeve expansion in certain cases. In this paper, Independent Component Analysis is exploited to linearly transform the random variables used in Karhunen–Loeve expansion resulting into a set of random variables exhibiting higher order decorrelation. The stochastic wave equation is investigated for numerical illustration whereby the random stiffness coefficient is modeled as a non-Gaussian stochastic process. Under the assumption of independence among the random variables used in the Karhunen–Loeve expansion and Independent Component Analysis representations, the latter provides more accurate statistical characterization of the output process for the specific cases examined. & 2014 Elsevier Ltd. All rights reserved.

1. Introduction Karhunen–Loeve Expansion (KLE) is the mean-square optimal representation of a random process using the spectral decomposition (e.g. [1,2]). KLE is most effective in capturing dominant modes of fluctuations of a random process by a few uncorrelated random variables (e.g. [1,2]). However, these random variables are usually not statistically independent for non-Gaussian processes. Independent Component Analysis (ICA) is a statistical estimation technique which provides a representation of a set of statistically dependent random variables by a linear combination of independent random variables,

Abbreviations: ICA, Independent Component Analysis; KLE, Karhunen–Loeve expansion; MLE, Maximum Likelihood Estimation; PDE, partial differential equation; SPDE, stochastic partial differential equation; POD, Proper Orthogonal Decomposition; pdf, probability density function; cdf, cumulative distribution function; MPI, Message Passing Interface; MCS, Monte Carlo Simulation n Corresponding author. Tel.: þ1 613 520 2600x6320; fax: þ 1 613 520 3951. E-mail address: [email protected] (A. Sarkar). 1 Currently at Sandia National Laboratories, P.O. Box 969 MS 9051, Livermore, CA 94551-9051, United States. http://dx.doi.org/10.1016/j.jsv.2014.04.027 0022-460X/& 2014 Elsevier Ltd. All rights reserved.

M. Khalil, A. Sarkar / Journal of Sound and Vibration 333 (2014) 5600–5613

5601

namely independent components [3–7]. Using a few realizations of the underlying dependent random variables, ICA identifies a linear transformation matrix that maps the original random variables to new components with higher order decorrelation properties. In the context of a stochastic partial differential equation (SPDE) with random coefficients, the random scatter of the system coefficients may be modeled as stochastic processes whenever substantial data is available. For computational efficiency, these random processes may be represented using KLE in terms of a few random variables. Note that KLE has been used in the field of structural dynamics for model reduction and uncertainty quantification (e.g. [8–20]). The output of the system may therefore be efficiently estimated using various uncertainty quantification techniques (e.g. [2]). The dependence among the KLE random variables may be ignored in many practical cases without incurring significant errors in the output statistics. This simplifying assumption may lead to significant inaccuracies in the statistics of the system response in some cases. In this investigation, ICA will be used to provide an improved representation of stochastic coefficients of the system where the representative random variables are as independent as possible. Of course, this representation will inherit the efficiency of KLE but reduces the error when the assumption of independence among the random variables is invoked. A similar concept has been reported for dynamic electric circuit simulation [21]. In this paper, our initiative is directed towards the application of ICA for uncertainty quantification in stochastic computational mechanics. The preliminary results on applications of ICA for non-Gaussian process representation are reported by the authors in the following conferences [22,23]. In this paper, the accuracy of input representation using ICA is investigated on the output of a SPDE describing wave propagation through random media. The random scatter in the stiffness coefficient of the media is modeled as a non-Gaussian random process. An analytical KLE representation of the stiffness coefficient is used and the associated output statistics of the system is estimated using Monte Carlo Simulation (MCS). The implication of the assumption of independence among the KLE random variables on the output statistics is then investigated. ICA is then applied to obtain an improved representation of the stiffness process whereby the random variables possess higher order decorrelation properties. Consequently, the random variables in the ICA representation are simulated under the assumption of statistical independence for computational efficiency. These three simulation results are compared. In this paper, KLE and ICA representations of uncertain system parameters are presented in Section 2. Next, the Maximum Likelihood Estimation (MLE) algorithm for ICA is reviewed in Section 3. Section 4 reports the results of numerical investigations which compares the accuracy of KLE and ICA representations applied to stochastic systems. In Section 5, the paper summarizes the findings of the paper. 2. Mathematical representation of stochastic processes Consider the following stochastic partial differential equation (PDE) that governs the behavior of numerous engineering systems: Lðuðx; θÞ; αðx; θÞÞ ¼ f ðx; θÞ;

(1)

where L, u, α and f are the random differential operator, solution process, random system parameter and deterministic input respectively, x and θ define the physical coordinate and random events in the probability space respectively. The approximation of a zero mean process α by KLE is given by (e.g. [1,2]) L

αðx; θÞ ¼ αðxÞ þ ∑ ξj ðθÞφj ðxÞ

(2)

j¼1

pffiffiffiffi L where αðxÞ denotes the mean, φj ¼ λj f j ðxÞ and fλj ; f j gj ¼ 1 are the first L-th eigenvalues and eigenvectors of the covariance function Rαα ðx; yÞ of αðx; θÞ respectively and the number of terms L is optimally selected to capture most (e.g. 99%) of the energy of the process; and fξj gLj ¼ 1 are uncorrelated, but not necessarily independent random variables whose joint probability density function (pdf) may be difficult to obtain in many cases. In this paper, we exploit ICA [22,24,25] in order to identify a set of random variables fsj gLj ¼ 1 having higher-order decorrelation property (in contrast to fξj gLj ¼ 1 ) using the following linear transformation: s ¼ Bξ:

(3)

In this investigation, B is obtained by the MLE algorithm outlined in [24,25] and described in the next section. Using ICA, the approximation of the stochastic process α becomes L

αðx; θÞ ¼ αðxÞ þ ∑ sj ðθÞψ j ðxÞ:

(4)

j¼1

The so-called ICA modes ψ j ðxÞ are given by L

ψ j ðxÞ ¼ ∑ aji φi ðxÞ;

(5)

i¼1

where the elements of the inverse of the matrix B define aji. Due to higher order decorrelation properties of the elements of fsg, the product of the marginal pdfs provides a reasonable approximation of the joint pdf of fsg.

5602

M. Khalil, A. Sarkar / Journal of Sound and Vibration 333 (2014) 5600–5613

In the context of structural reliability analysis, the Rosenblatt transformation is generally used to map non-Gaussian random vectors into a set of standard normal random variables through nonlinear transformations (e.g. [26]). The usefulness of the Rosenblatt transformation in tackling practical applications is generally limited as it requires the joint or conditional probability density functions of the non-Gaussian random vectors [26]. When the available data is sparse, the joint or conditional probability density functions of a non-Gaussian random vector are difficult to estimate with sufficient accuracy. In contrast to the Rosenblatt transformation, the transformation based on ICA does not require any prior knowledge of the joint or conditional probability density functions of the (non-Gaussian) KLE random variables. The probabilistic characterizations of ICA random variables are obtained directly from data by satisfying higher order decorrelation properties without the prior knowledge of the joint or conditional probability densities of the KLE random variables.

3. ICA by MLE Closely following the books by Roberts and Everson [24] and Hyvärinen et al. [25] and using the well-known results of linear transformation [27], the maximum likelihood estimation procedure of the mixing matrix B is reviewed for completeness. For the linear mapping ξ ¼ B  1 s;

(6)

the joint pdf of ξ can be written as [24,25] pðξÞ ¼ jdet BjpðsÞ L

¼ jdet Bj ∏ pi ðsi Þ i¼1 L

T

¼ jdet Bj ∏ pi ðbi ξÞ;

(7)

i¼1

where pi and bi are the pdf of the independent component si and the i-th row of B respectively. In the current investigation, we make the following assumptions which ensure that the ICA model is identifiable (and thus the mixing matrix B is fully ranked) [24,25]: (a) the independent components are statistically independent, (b) the independent components are non-Gaussian, and (c) the number of independent components is the same as that of observed variables (KLE modes), leading to a square mixing matrix. In maximum likelihood estimation, two approximations of the densities pi of an independent component provide locally consistent estimates [24,25] (see Appendix A). Assuming that N observations of ξ, denoted by ξð1Þ; …; ξðNÞ, are available, the likelihood function can be obtained as [24,25] N

L

k¼1

i¼1

T

LðBÞ ¼ ∏ jdet Bj ∏ pi ðbi ξðkÞÞ:

(8)

In Appendix B, the gradient of the log-likelihood function with respect to the matrix B is shown to be [24,25]   1 ∂ ln L ¼ ðBT Þ  1 þE gðBξÞξT N ∂B

(9)

where EðÞ denotes the expectation operator, gi, the i-th element of g, defines the score function of the pdf of si. To accelerate the convergence of MLE algorithms based on iterative methods, the above expression of the gradient can be modified using the natural gradient as follows [24,25]:   B’B þ μðI þ E gðsÞξT BT ÞB; (10) where μ is the step length which is obtained optimally by performing one-dimensional line search optimization using the golden section bracketing  method  [28].  The condition E gðsÞξT BT ¼ E gðsÞsT ¼ I dictates the convergence of the iterative algorithm which can be interpreted as the higher order (nonlinear) decorrelation of fsg and fgðsÞg. ICA has some difficulties and relies on several assumptions [24,25]. The obvious but non-trivial difficulty is encountered when estimating a nonlinear mixing model. One inherent difficulty in nonlinear ICA relates to the existence and uniqueness of the solution if the space of nonlinear mixing models is too small or too large respectively [24,25]. Another difficulty stems from the computational demands required by the available methods to the nonlinear ICA problem [24,25]. The ICA model discussed in this investigation is a linear model, in that the dependent random variables are obtained by a linear transformation of independent random variables. There are several proposed algorithms that strive to estimate a nonlinear ICA model. For example, one can use an extension to the previously presented ICA model to include certain forms of nonlinearities (e.g. [29]). Another difficulty is faced when estimating under-determined or over-determined models. Such models occur when the number of independent random variables is not equal to that of dependent random variables.

M. Khalil, A. Sarkar / Journal of Sound and Vibration 333 (2014) 5600–5613

5603

4. Application to stochastic partial differential equations The following partial differential equation describes one dimensional axial wave propagation in a heterogeneous medium (e.g. [13,14,30–32],)   ∂ ∂u ∂2 u ∂2 u ∂u EðxÞAðxÞ þ γ 1 ðxÞ ¼ ρðxÞAðxÞ 2 þγ 2 ðxÞ (11) ∂x ∂x ∂t∂x ∂t ∂t where E is the modulus of elasticity, ρ is the density of the medium, A is the cross sectional area, γ1 and γ2 are the damping coefficients per unit length and u represents the dynamic motion. The random heterogeneity in the system properties (e.g. the modulus of elasticity, cross sectional area, damping coefficients) will render Eq. (11) into a time dependent stochastic partial differential equations (e.g. [13,14,31,32]). The finite element discretization of a one-dimensional second order (with respect to spatial derivatives) partial differential equation leads to a coupled array of mass–spring oscillators (e.g. [33]) as shown in Fig. 1. For an n degree-of-freedom model obtained from discretization of SPDE, the mass and stiffness matrices of the system are given by 2 3 m 0 60 m ⋱ 7 6 7 M¼6 (12) 7; 4 ⋱ ⋱ 05 0 2

k1 þ k2 6 6  k2 K¼6 6 4

m 3

k2 k2 þk3

⋱

⋱

⋱

 kn

 kn

kn þ kn þ 1

7 7 7; 7 5

(13)

and the damping matrix is assumed to be C ¼ α0 M þα1 K where α0 and α1 are the proportional damping coefficients. The system is fixed at both ends.

Fig. 1. Linear array of mass–spring oscillators (adapted from [15]).

mean impulse response

x 10−8 1 0 −1 0

0.02

0.04

0.06

0.08

0.1 t

0.12

0.14

0.16

0.18

0.02

0.04

0.06

0.08

0.1 t

0.12

0.14

0.16

0.18

−9

standard deviation

x 10 6 4 2 0 0

Fig. 2. Driving-point impulse response function of a coupled linear array of mass–spring oscillators at node 50 with stiffness modeled as a truncated Taylor series expansion of a log-normal process: (a) mean response, (b) standard deviation.

5604

M. Khalil, A. Sarkar / Journal of Sound and Vibration 333 (2014) 5600–5613

For simple analytical representation, the stiffness coefficients are modeled as a truncated Taylor series expansion of a lognormal process. For this process, one can find an exact KLE representation as described later. The random stiffness coefficients are given by ki ¼ k½1 þ ϵαðxi ; θÞ;

(14)

where k is the mean stiffness, ϵ is the coefficient of variation and αðxi ; θÞ is a truncated Taylor series expansion of a lognormal process (i.e. expðgðx; θÞÞ having zero mean and unit variance expressed as qffiffiffiffiffiffiffiffiffi 320  5=8 þ g ðxi ; θÞ þg 2 ðxi ; θÞ=2 αðxi ; θÞ ¼ 1341 þ g 3 ðxi ; θÞ=6 þ g 4 ðxi ; θÞ=24 þg 5 ðxi ; θÞ=120;

(15)

with gðxi ; θÞ ¼ AðθÞ cos ðωxi Þ þ BðθÞ sin ðωxi Þ:

(16)

A and B define independent Gaussian random variables with zero mean and unit standard deviation. The truncated KLE representation of αðxi ; θÞ is ! rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi" 1600 2A2 2B2 A4 A2 B2 B4 5 þ þ þ þ  αðxi ; θÞ ¼ 199 809 3 3 24 12 24 3 ! 3 2 5 3 2 4 8A A AB A A B AB þ þ þ þ þ þ cos ðωxi Þ 3 3 3 72 36 72 ! 8B B3 A2 B B5 A2 B3 A4 B þ þ þ þ þ þ sin ðωxi Þ 3 3 72 3 72 36 ! 2A2 2B2 A4 B4  þ  þ cos ð2ωxi Þ 3 3 18 18 ! 4AB A3 B AB3 þ þ þ sin ð2ωxi Þ 3 9 9 ! A3 AB2 A5 A3 B2 AB4  þ   þ cos ð3ωxi Þ 9 3 144 72 48 ! B3 A2 B B5 A2 B3 A4 B  þ þ sin ð3ωxi Þ þ  þ 3 48 9 144 72 ! A4 A2 B2 B4  þ þ cos ð4ωxi Þ 72 12 72

−8

mean impulse response

x 10 1 0 −1

0

0.02

0.04

0.06

0.08

0.1 t

0.12

0.14

0.16

0.18

0.02

0.04

0.06

0.08

0.1 t

0.12

0.14

0.16

0.18

−9

standard deviation

x 10 4 2 0 0

Fig. 3. Cross impulse response function of a coupled linear array of mass–spring oscillators between nodes 25 and 75 with stiffness modeled as a truncated Taylor series expansion of a log-normal process: (a) mean response, (b) standard deviation.

M. Khalil, A. Sarkar / Journal of Sound and Vibration 333 (2014) 5600–5613

þ

A3 B AB3  18 18

5605

! sin ð4ωxi Þ

! A5 A3 B2 AB4  þ cos ð5ωxi Þ 720 72 144 # ! B5 A2 B3 A4 B  þ þ sin ð5ωxi Þ : 144 720 72 þ

(17)

The KLE random coefficients given by the amplitudes of the harmonics are statistically uncorrelated, but dependent nonGaussian random variables. As an example, the amplitudes of the cos ð4ωxi Þ and sin ð4ωxi Þ terms are uncorrelated since " ! !# h  i A4 A2 B2 B4 A3 B AB3 1  E A7 B  AB7 þ 7 A5 B3  A3 B5  þ E ¼ 18 1296 72 12 72 18 ¼ 0; but they are statistically dependent being the functions of the same random variables A and B. Based on the KLE representation in Eq. (17), the exact mixing model is given by the nonlinear mapping: ! 9 8 rffiffiffiffiffiffiffi 9 2A2 2B2 A4 A2 B2 B4 5 > > > > > > þ þ þ þ  > > > > 37 3 3 3 24 12 24 > > > > > > rffiffiffiffiffiffiffiffiffiffi > > ! > > 3 2 5 3 4 > > 2 > > 27 8A A AB A A B AB > > > > þ þ þ þ þ > > > > > > 562 3 3 3 72 36 72 > > > > > > ! > > r ffiffiffiffiffiffiffiffiffi ffi > > 2 2 3 4 3 5 > > > 27 8B B A B B A B A B > > > > > þ þ þ þ þ > > > 562 3 3 72 > 3 72 36 > > > > > > > > ! ffiffiffiffiffiffiffiffiffi ffi r > > > > 2 4 2 4 > > 27 2A 2B A B > > > > > >  þ  > > > > 112 3 3 18 18 > > > > > > ! > > r ffiffiffiffiffiffiffiffiffi ffi > > 3 3 > > > > 27 4AB A B AB > > > > þ þ > > > > > > 112 3 9 9 > > > > > > ! > > r ffiffiffiffiffiffi ffi > >

 3 2 5 3 4 2 = < 54 A AB A A B AB A  þ   : ⟼ 37 9 3 144 72 48 > > B > > > > > rffiffiffiffiffiffiffi ! > > > > > > > 54 B3 A2 B B5 A2 B3 A4 B > > > >  þ þ  þ > > > > > > 37 3 48 9 144 72 > > > > > > ! > > > > 4 2 2 4 > > p ffiffiffiffiffiffi > > A A B B > > > >  þ 27 > > > > 72 12 72 > > > > > > > > ! > > > > 3 3 > > pffiffiffiffiffiffi A B AB > > > > > >  27 > > > > 18 18 > > > > > > ! > > > > 5 3 4 2 > > pffiffiffiffiffiffiffiffiffi A > > A B AB > > > >  þ 270 > > > > > > 720 72 144 > > > > > > ! > > > > 2 3 4 5 > > p ffiffiffiffiffiffiffiffi ffi > > B A B A B > > > >  þ 270 > > ; : 144 720 72

(18)

(19)

Note that the underlying mixing model is nonlinear and over-determined as the dimension of the dependent random vector exceeds that of the independent counterpart. This nonlinear model is conceptually and computationally difficult to estimate [24,25]. In the linear ICA formalism, the dimension of the dependent and independent components are the same. This model is identifiable. Although the components obtained may not be dependent, they will satisfy higher order decorrelation properties. The following values are used in the numerical experiments: n¼100, α0 ¼ 14, α1 ¼ 0, m¼ 0.1, k ¼ 3:67  106 , ϵ ¼ 0:1, ω ¼ 2π and xi ¼ ðð2i 1Þ=2ðn þ 1ÞÞL for i ¼ 1; 2; …; n þ 1. where L denotes the length of the domain over which the SPDE in Eq. (11) is defined. The numbers of finite elements and nodes are taken to be 101 and 102 respectively, leading to the total degree-of-freedom n ¼100. Note that xi denotes the location of the mid-point of the i-th element. For simplicity, the element stiffness matrix is calculated analytically using the value of the stochastic process αðxi ; θÞ at the mid-point of the element. For these numerical values of the mean k and the coefficient of variation ϵ, the probability of a negative stiffness coefficient sample is negligibly small. The mean driving-point impulse response function of node 50 and the cross impulse response function between nodes 75 and 25, along with the associated standard deviations, are shown in Figs. 2 and 3 respectively. Although other statistical quantities (e.g. higher order moments or entropy) are more appropriate measure of uncertainty for non-Gaussian output such as the impulse response function, the standard deviation is considered to be

5606

M. Khalil, A. Sarkar / Journal of Sound and Vibration 333 (2014) 5600–5613

adequate to highlight the usefulness of the current investigation. These results are used as a benchmark to test the accuracy of both KLE and ICA representations. Next, the stiffness process is approximated by its KLE expansion given by Eq. (17) ignoring the statistical dependence among the amplitude terms. Using ICA, the random vector, whose elements are these random amplitude terms, is transformed using MLE of the inverse mixing B obtained from the algorithm in Eq. (10) as 2 3  0:212  0:229 0:291 0:045 0:398 0:396 0:244 0:480  0:090 0:365  0:261 6 0:348 0:469 0:171 0:369 0:320 0:163 0:391 0:035 0:362 0:204 0:195 7 6 7 6 7 6  0:358  0:146 7 0:044 0:336 0:281 0:196 0:190  0:234 0:227  0:312 0:615 6 7 6  0:222 0:105 0:401 0:385 0:127  0:084  0:365 0:244 0:218 0:566 0:216 7 6 7 6 7 6 0:303 0:164 0:445  0:302 0:365  0:449  0:010 0:247  0:334 0:112  0:268 7 6 7 6 7  0:265 0:348  0:177 0:409 0:413 0:134 0:382 0:274 0:053  0:364 7: (20) B ¼ 6 0:254 6 7 6  0:316 0:234 0:306 0:098 0:398  0:422 0:189 0:392 0:180 0:249  0:348 7 6 7 6 7  0:437 0:148 0:352 0:282  0:278  0:347 0:107 0:391 0:111  0:356 7 6 0:287 6 7 6 0:350  0:439 0:170 0:370 0:289  0:217 0:342 0:123  0:381 0:045 0:330 7 6 7 6 7 0:400 0:115 0:330 0:242 0:293  0:340 0:181  0:422 0:063  0:364 5 4 0:333  0:289 0:009 0:368 0:363 0:057  0:112 0:440 0:487  0:145 0:173  0:389

0.5

0

−2

−1

0 1 α( x10)

0.5

0

2

−2

−1

0

1

0

2

0 1 α( x )

−1

1

2

0 1 α( x )

2

1.2

1

60

0.8 0.6

1 0.8 0.6

0.4

0.4

0.2

0.2 −2

2

0 30

p(α( x ))

50

p(α( x ))

p(α( x )) 40

0.5

−1

−2

α( x )

1.2

−2

0.5

α( x20)

1

0

1

30

p(α( x ))

1

20

p(α( x ))

1

10

p(α( x ))

Note that the matrix B is an orthogonal matrix. As the KLE random coefficients are decorrelated, the ICA random coefficients are also decorrelated. Furthermore, the higher order correlations among the random variables in the ICA

−1

0

1

−2

2

−1

α( x ) 50

40

60

1.4 1

80

0.8 0.6

1

90

p(α( x ))

1

p(α( x ))

p(α( x70))

1.2

0.5

0.5

0.4 0.2 −2

−1

0 α( x70)

1

2

0

−2

−1

0 1 α( x80)

2

0

−2

−1

0

1

2

α( x90)

Fig. 4. pdfs of the stiffness process using the original representation and KLE and ICA approximate representations assuming independence among their respective random coefficients. The ICA model was obtained using 50 samples of the KLE random coefficients. The stiffness process is a truncated Taylor series expansion of a log-normal process (solid line), and the KLE (solid line with stars) and ICA (solid line with circles) approximations of stiffness process assuming independence among the random coefficients.

M. Khalil, A. Sarkar / Journal of Sound and Vibration 333 (2014) 5600–5613

5607

0.5

0

−2

−1

0 1 α( x10)

0.5

0

2

1

30

p(α( x ))

1

20

p(α( x ))

1

10

p(α( x ))

representation are smaller in relation to the KLE random variables. Therefore, the stiffness process can thus be approximated efficiently and accurately by its ICA expansion ignoring the statistical dependency among the ICA random coefficients. The pdfs of the original stiffness coefficient ki for i¼ 20 and its KLE and ICA representations are shown in Figs. 4 and 5 whereby the ICA MLE estimate of the inverse mixing matrix B was obtained using 50 and 1000 random samples of the KLE coefficients in Eq. (17), respectively. The KLE and ICA random variables are sampled while ignoring the statistical dependence among the random variables. To obtain accurate estimates of the pdfs of the stiffness coefficient, 5  106 Monte Carlo samples are used for each representation. Such a large number of samples were required to obtain accurate pdf estimates (but not needed to estimate B). Under the assumption of independence of random variables, the ICA representation of the stiffness process provides more accurate pdfs of the stiffness coefficients in contrast to the KLE representation. When 1000 samples of the KLE random variables are used to estimate the ICA inverse mixing matrix B (see Fig. 5) as opposed to 50 random samples (see Fig. 4), the ICA representation provides more accurate marginal pdfs of the stiffness coefficients. For the subsequent results, we will use the ICA representation that is obtained using 1000 random samples of the KLE coefficients. The seven KLE modes φj and ICA modes ψj of the process αðxi ; θÞ are shown in Fig. 6. The MLE estimate of the ICA inverse mixing matrix B is obtained using 1000 samples of the KLE random coefficients. The KLE modes φj are the harmonics in Eq. (17). As noted earlier, based on the linear ICA model assumed in Eq. (3), the ICA modes are the linear superposition of KLE modes. The coefficients aij mapping the KLE modes to the ICA modes in Eq. (5) are plotted in Fig. 6. A measure of uncertainty in the estimates of these coefficients is not directly provided by the MLE algorithm. The mean of the driving-point and cross impulse response functions using the original stiffness process in Eq. (14), the KLE and ICA approximations are shown in Figs. 7a and 8a, respectively. For the KLE and ICA representations, the random

−2

−1

0 1 α( x20)

0.5

0

2

−2

−1

0

1

2

0 1 α( x )

2

0 1 α( x90)

2

α( x30)

1.4 1.2

0.5

0

−2

−1

0 1 α( x )

1

60

p(α( x ))

50

p(α( x ))

p(α( x )) 40

1.2 1

0.8 0.6

1 0.8 0.6

0.4

0.4

0.2

0.2 −2

2

−1

0 1 α( x )

−2

2

−1

50

40

60

1

80

0.8 0.6 0.4

1

90

p(α( x ))

1

p(α( x ))

p(α( x70))

1.2

0.5

0.5

0.2 −2

−1

0 1 α( x70)

2

0

−2

−1

0 1 α( x80)

2

0

−2

−1

Fig. 5. pdfs of the stiffness process using the original representation and KLE and ICA approximate representations assuming independence among their respective random coefficients. The ICA model was obtained using 1000 samples of the KLE random coefficients. The stiffness process is a truncated Taylor series expansion of a log-normal process (solid line), and the KLE (solid line with stars) and ICA (solid line with circles) approximations of stiffness process assuming independence among the random coefficients.

5608

M. Khalil, A. Sarkar / Journal of Sound and Vibration 333 (2014) 5600–5613

0.5 0 −0.5

1 0 −1 0

1

2

3 0.5 0 −0.5

1 0 −1 0

1

2

3 0.5 0 −0.5

1 0 −1 0

1

2

3 0.5 0 −0.5

1 0 −1 0

1

2

3 0.5 0 −0.5

1 0 −1 0

1

2

3 0.5 0 −0.5

1 0 −1 0

1

2

3 0.5 0 −0.5

1 0 −1 0

1

2

3

1 0 −1 a a a a a a a a a a a

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

1 0 −1 a a a a a a a a a a a

1 0 −1 a a a a a a a a a a a

1 0 −1 a a a a a a a a a a a

1 0 −1 a a a a a a a a a a a

1 0 −1 a a a a a a a a a a a

1 0 −1 a a a a a a a a a a a

Fig. 6. (a–g) KLE modes φj for j ¼ 1; …; 7 of the stochastic stiffness process; (h–n) coefficients aij of the maximum likelihood estimate of the mixing matrix relating the KLE modes to the ICA modes; (o–u) ICA modes ψi for i ¼ 1; …; 7 of the stochastic stiffness process. Note that the KLE and ICA modes are plotted for x between 0 and 3.

coefficients ξi and si are sampled independently ignoring their statistical correlation which offers computational efficiency. The standard deviations of the impulse response functions are shown in Figs. 7b and 8b using 15,000 realizations of the stiffness process executed on an Intel Xeon cluster with InfiniBand interconnect using Message Passing Interface (MPI) [34]. The ICA MLE algorithm is applied using 1000 samples of KLE coefficients to estimate the inverse mixing matrix B. Both KLE and ICA representations provide reasonably good estimates of the mean impulse response functions when the dependence among the random variables is ignored with ICA showing slightly greater accuracy. However, the ICA representation provides improved estimates of the standard deviation in contrast to the KLE representation. The assumption of independence of random variables in the KLE representation most severely affects the standard deviation estimates of the cross impulse response function, leading to normalized errors of approximately 50% in some cases in comparison to the baseline case. An accurate characterization of cross impulse response function depends critically on the reliable representation of the randomness in the stiffness coefficient. The pdfs of the driving-point and cross impulse response functions using the original stiffness process, and KLE and ICA approximations at t ¼0.0833 seconds are shown in Figs. 9 and 10, respectively. The associated cumulative distribution functions (cdfs) are shown in Figs. 11 and 12. This specific time was chosen since the standard deviation of the impulse response functions is at its peak. The strong non-Gaussian characteristics are evident. Clearly both the KLE and ICA representations (with the assumption of independence among the coefficients) lead to approximation errors. However ICA

M. Khalil, A. Sarkar / Journal of Sound and Vibration 333 (2014) 5600–5613

5609

−8

mean impulse response

x 10 1 0

Original KLE ICA

−1 0

0.02

0.04

0.06

0.08

0.1 t

0.12

0.14

0.16

0.18

−9

standard deviation

x 10 6 4

Original KLE ICA

2 0 0

0.02

0.04

0.06

0.08

0.1 t

0.12

0.14

0.16

0.18

Fig. 7. Driving-point impulse response function of a coupled linear array of mass–spring oscillators at node 50 with the stiffness given by a truncated Taylor series expansion of a log-normal process (solid line), and the KLE (solid line with stars) and ICA (solid line with circles) approximations of stiffness process assuming independence among the random coefficients: (a) mean response and (b) standard deviation.

mean impulse response

x 10

8

1 0 Original KLE ICA

1 0

standard deviation

x 10

0.02

0.04

0.06

0.08

0.1 t

0.12

0.14

0.16

0.18

9

6 4 Original KLE ICA

2 0 0

0.02

0.04

0.06

0.08

0.1 t

0.12

0.14

0.16

0.18

Fig. 8. Cross impulse response function of a coupled linear array of mass–spring oscillators at node 50 with the stiffness given by a truncated Taylor series expansion of a log-normal process (solid line), and the KLE (solid line with stars) and ICA (solid line with circles) approximations of stiffness process assuming independence among the random coefficients: (a) mean response and (b) standard deviation.

outperforms KLE in this case in terms of maximum absolute error and root-mean-squared error in the representation of the pdf (although these results are not shown in the paper). 5. Conclusion When the assumption of independence is invoked, ICA representation of a non-Gaussian input process is superior to that of KLE. Hence ICA based approximate representation of non-Gaussian input parameters naturally provides accurate statistics of output processes. In this investigation, the non-Gaussian stiffness coefficient of the one-dimensional stochastic wave

5610

M. Khalil, A. Sarkar / Journal of Sound and Vibration 333 (2014) 5600–5613

x 107 Original KLE ICA

18 16 14

pdf

12 10 8 6 4 2 0

−1.5

−1

−0.5

0

0.5

1

1.5

2 −8

x 10

Impulse response

Fig. 9. PDF of driving-point impulse response function at 0.0833 s of a coupled linear array of mass–spring oscillators at node 50 with the stiffness given by a truncated Taylor series expansion of a log-normal process (solid line), and the KLE (solid line with stars) and ICA (solid line with circles) approximations of stiffness process assuming independence among the random coefficients.

x 10

8

Original KLE ICA

2

pdf

1.5

1

0.5

0

−1.5

−1

−0.5

0

0.5

1

1.5

Impulse response

2 −8

x 10

Fig. 10. PDF of cross impulse response function at 0.0833 s of a coupled linear array of mass–spring oscillators at node 50 with the stiffness given by a truncated Taylor series expansion of a log-normal process (solid line), and the KLE (solid line with stars) and ICA (solid line with circles) approximations of stiffness process assuming independence among the random coefficients.

equation is represented analytically by KLE. The random variables in KLE are statistically orthogonal (linearly decorrelated), but not independent. Next the non-Gaussian random variables of KLE are transformed to new variables which have higher order (nonlinear) decorrelation properties. The effect of the assumption of independence among the KLE and ICA random variables is then investigated on the output statistics (i.e. mean and standard deviation). Based on the numerical experiments, the following observations are noted:

 Under the assumption of independence of the random variables, ICA provides a more accurate representation of the stiffness process in contrast to KLE.

 Both ICA and KLE representations of the stiffness process lead to accurate mean estimates of the system output in relation to the baseline (exact) case, albeit with slightly greater accuracy for ICA.

 In contrast to KLE, ICA representation of the stiffness process provides more accurate standard deviation estimates of the system output under the assumption of independence among its random variables.

M. Khalil, A. Sarkar / Journal of Sound and Vibration 333 (2014) 5600–5613

1

5611

Original KLE ICA

0.8

cdf

0.6

0.4

0.2

0 −1.5

−1

−0.5

0

0.5

1

1.5

2 −8

x 10

Impulse response

Fig. 11. CDF of driving-point impulse response function at 0.0833 s of a coupled linear array of mass–spring oscillators at node 50 with the stiffness given by a truncated Taylor series expansion of a log-normal process (solid line), and the KLE (solid line with stars) and ICA (solid line with circles) approximations of stiffness process assuming independence among the random coefficients.

1

Original KLE ICA

0.8

cdf

0.6

0.4

0.2

0 −1.5

−1

−0.5

0 Impulse response

0.5

1

1.5

2 x 10−8

Fig. 12. CDF of cross impulse response function at 0.0833 s of a coupled linear array of mass–spring oscillators at node 50 with the stiffness given by a truncated Taylor series expansion of a log-normal process (solid line), and the KLE (solid line with stars) and ICA (solid line with circles) approximations of stiffness process assuming independence among the random coefficients.

The methodology reported in this paper can be exploited for an efficient spectral stochastic finite element [2] formulation involving polynomial chaos representation of non-Gaussian input process. The application of ICA in conjunction with polynomial chaos expansions for stochastic dynamical systems will be the subject of future investigation. Appendices The books by Roberts and Everson [24] and by Hyvärinen et al. [25] are followed closely in these Appendices. Appendix A. Densities of the independent components The likelihood function for estimating B involves pdf pi ðsi Þ of the independent component si which is a priori unknown. b iþ and subGaussian pdf p b i approximations of pi ðsi Þ However, it has been shown that the following two superGaussian pdf p

5612

M. Khalil, A. Sarkar / Journal of Sound and Vibration 333 (2014) 5600–5613

are sufficient to obtain consistent estimates of B [25,35] þ

b i ðsi Þ ¼ α1  2 ln coshðsi Þ ln p 

b i ðsi Þ ¼ α2  s2i =2 þ ln coshðsi Þ ln p

(A.1) (A.2)

where α1 and α2 are positive parameters. Appendix B. Gradient of the log-likelihood function Consider the log-likelihood function for B

( ) n 1 T ln LðBÞ ¼ lnjdet Bjþ E ∑ ln pi ðbi ξÞ N i¼1

(B.1)

for which the derivative of the second term in the right hand side with respect to the element brq can be written as [25] ( ) ( )  n ∂ ∂ n T T E ∑ ln pi ðbi ξÞ ¼ E ∑ ln pi bi ξ ∂brq i ¼ 1 ∂brq i ¼ 1 ( !) n ∂ ∑ ln pi ∑bij ξj ¼E ∂brq i ¼ 1 j ( !) ∂ ¼E ln pr ∑brj ξj ∂brq j

 ∂ ln pr ðsr Þ ¼E ∂brq

 ∂ ln pr ðsr Þ ∂sr ¼E ∂sr ∂brq ( ) T ∂b ξ ¼ E g r ðsr Þ r ∂brq 9 8 brj ξj > > = <  ∂∑ j T ¼ E g r br ξ > ∂b rq > ; : n o T ¼ E g r ðbr ξÞξq : (B.2) In matrix form, we obtain (

∂ E ∂B

)

n

T

∑ ln pi ðbi ξÞ

i¼1

9 9 88 > > > g 1 ðbT1 ξÞ > > > > = = <> < ⋮ ½ξ1 ⋯ ξn  ¼E > > > > > > ; ; :> : g n ðbT ξÞ > n

¼ EfgðBξÞξT g:

(B.3)

where ∂ b i ðsi Þ ln p ∂si b i ðsi Þ 1 ∂p : ¼ b i ðsi Þ ∂si p

g i ðsi Þ ¼

Using the identity

∂ ∂B

(B.4)

  ln det BjÞ ¼ ðBT Þ  1 , the gradient of the log-likelihood function can be expressed as [24,25]   1 ∂ ln L ¼ ðBT Þ  1 þE gðBξÞξT : N ∂B

(B.5)

References [1] P. Holmes, J.L. Lumley, G. Berkooz, Turbulence, Coherent Structures, Dynamical Systems and Symmetry, Cambridge University Press, New York, 1996. [2] R.G. Ghanem, P.D. Spanos, Stochastic Finite Elements: A Spectral Approach, Springer-Verlag, 1991. [3] M. Hämäläinen, R. Hari, R. Ilmoniemi, J. Knuutila, O.V. Lounasmaa, Magnetoencephalography—theory, instrumentation, and applications to noninvasive studies of the working human brain, Reviews of Modern Physics 65 (2) (1992) 413–497. [4] K. Kiviluoto, E. Oja, Independent Component Analysis for parallel financial time series, Proceedings of the ICONIP'98, International Conference on Neural Information Processing, Vol. 2, Tokyo, Japan, 1998, pp. 895–898. [5] A. Hyvärinen, P.O. Hoyer, M. Inki, Topographic independent component analysis, Neural Computation 13 (2001).

M. Khalil, A. Sarkar / Journal of Sound and Vibration 333 (2014) 5600–5613

5613

[6] Z.N. Li, Y.Y. He, F.L. Chu, J. Han, W. Hao, Fault recognition method for speed-up and speed-down process of rotating machinery based on independent component analysis and factorial hidden Markov model, Journal of Sound and Vibration 291 (1–2) (2006) 60–71. [7] M.J. Zuo, J. Lin, X.F. Fan, Feature separation using ICA for a one-dimensional time series and its application in fault detection, Journal of Sound and Vibration 287 (3) (2005) 614–624. [8] Philipp Glüsmann, Edwin Kreuzer, On the application of Karhunen–Loeve transform to transient dynamic systems, Journal of Sound and Vibration 328 (4–5) (2009) 507–519. [9] S. Adhikari, M.I. Friswell, Distributed parameter model updating using the Karhunen–Loeve expansion, Mechanical Systems and Signal Processing 24 (2) (2010) 326–339. [10] Xianghong Ma, Alexander F. Vakakis, Lawrence A. Bergman, Karhunen-Loeve analysis and order reduction of the transient dynamics of linear coupled oscillators with strongly nonlinear end attachments, Journal of Sound and Vibration 309 (3–5) (2008) 569–587. [11] Sergio Bellizzi, Rubens Sampaio, POMs analysis of randomly vibrating systems obtained from Karhunen–Loève expansion, Journal of Sound and Vibration 297 (3–5) (2006) 774–793. [12] A. Sarkar, R. Ghanem, Mid-frequency structural dynamics with parameter uncertainty, Computer Methods in Applied Mechanics and Engineering 191 (47) (2002) 5499–5513. [13] R. Ghanem, A. Sarkar, Reduced models for the medium-frequency dynamics of stochastic systems, Journal of the Acoustical Society of America 113 (2003) 834. [14] A. Sarkar, R. Ghanem, A substructure approach for the midfrequency vibration of stochastic systems, Journal of the Acoustical Society of America 113 (2003) 1922. [15] Mohammad Khalil, Sondipon Adhikari, Abhijit Sarkar, Linear system identification using proper orthogonal decomposition, Mechanical Systems and Signal Processing 21 (8) (2007) 3123–3145. [16] M. Amabili, A. Sarkar, M.P. Paıdoussis, Reduced-order models for nonlinear vibrations of cylindrical shells via the proper orthogonal decomposition method, Journal of Fluids and Structures 18 (2) (2003) 227–250. [17] M.P. Paıdoussis, A. Sarkar, C. Semler, A horizontal fluid-conveying cantilever: spatial coherent structures, beam modes and jumps in stability diagram, Journal of sound and vibration 280 (1) (2005) 141–157. [18] A. Sarkar, M.P. Paıdoussis, A compact limit-cycle oscillation model of a cantilever conveying fluid, Journal of Fluids and Structures 17 (4) (2003) 525–539. [19] A. Sarkar, M.P. Paidoussis, A cantilever conveying fluid: coherent modes versus beam modes, International Journal of Non-Linear Mechanics 39 (3) (2004) 467–481. [20] M. Amabili, A. Sarkar, M.P. Païdoussis, Chaotic vibrations of circular cylindrical shells: Galerkin versus reduced-order models via the proper orthogonal decomposition method, Journal of Sound and Vibration 290 (3) (2006) 736–762. [21] B. Wu, J. Zhu, F.N. Najm, Dynamic-range estimation, IEEE Transactions On Computer-Aided Design Of Integrated Circuits and Systems 26 (9) (2006) 1618–1636. [22] M. Khalil, A. Sarkar, Independent component analysis for uncertainty representation of stochastic systems, Proceedings of the 47th AIAA/ASME/ASCE/ AHS/ASC Structures, Structural Dynamics, and Materials Conference Schaumburg, USA, 2008. [23] M. Khalil, A. Sarkar, Karhunen–Loeve expansion versus independent component analysis for non-Gaussian stochastic processes: application to uncertain systems, Presented at 10th US National Congress on Computational Mechanics Conference, Ohio, USA, 2009. [24] S. Roberts, R. Everson, Independent Component Analysis: Principles and Practice, Cambridge University Press, Cambridge, 2001. [25] A. Hyvärinen, J. Karhunen, E. Oja, Independent Component Analysis, Wiley, New York, 2001. [26] Robert E. Melchers, Structural Reliability Analysis and Prediction, Civil Engineering Series, Wiley, Chichester, 1999. [27] W. Feller, An Introduction to Probability Theory and Its Applications, third edition, Wiley, New York, 1968. [28] J.A. Snyman, Practical Mathematical Optimization: An Introduction to Basic Optimization Theory and Classical and New Gradient-Based Algorithms, Springer, Heidelberg, 2005. [29] G. Burel, Digital holography and Karhunen–Loeve decomposition for the modal analysis of two-dimensional vibrating structures, Neural Network 5 (6) (1992) 937–947. [30] R.W. Clough, J. Penzien, Dynamics of structures, Vol. 634, McGraw-Hill, New York, 1993. [31] C.S. Manohar, S. Adhikari, Dynamic stiffness of randomly parametered beams, Probabilistic Engineering Mechanics 13 (1) (1998) 39–51. [32] C.S. Manohar, S. Adhikari, Statistics of vibration energy flow in randomly parametered trusses, Journal of Sound and Vibration 217 (1) (1998) 43–74. [33] G. Strang, Introduction to Applied Mathematics, Wellesley-Cambridge Press, Wellesley, MA, 1986. [34] Message Passing Interface 〈http://www.mpi-forum.org/〉. [35] Lei Xu, One-bit-matching ica theorem, convex-concave programming, and combinatorial optimization, in: Jun Wang, Xiaofeng Liao, Zhang Yi (Eds.), Advances in Neural Networks—ISNN 2005, Lecture Notes in Computer Science, Vol. 3496, Springer, Berlin, Heidelberg, 2005, pp. 5–20.