Orthogonality

The core idea of polynomial chaos expansions is that the polynomials used as an expansion are all mutually orthogonal. The relation is typically written mathematically as:

\[\left\langle \Phi_n, \Phi_m \right\rangle = 0 \qquad n \neq m\]

In practice this relation is instead expressed by the equivalent notation using expected values:

\[\mbox E\left(\Phi_n \Phi_m\right) = 0 \qquad n \neq m\]

In chaospy this property can be tested by taking the outer product of two expansions, and checking if the expected value of the resulting matrix is diagonal. For example, for a basic monomial:

import chaospy

monomial = chaospy.monomial(4)
monomial

polynomial([1, q0, q0**2, q0**3])

monomial2 = chaospy.outer(monomial, monomial)
monomial2

polynomial([[1, q0, q0**2, q0**3],
            [q0, q0**2, q0**3, q0**4],
            [q0**2, q0**3, q0**4, q0**5],
            [q0**3, q0**4, q0**5, q0**6]])

normal = chaospy.Normal(0, 1)
chaospy.E(monomial2, normal)

array([[ 1.,  0.,  1.,  0.],
       [ 0.,  1.,  0.,  3.],
       [ 1.,  0.,  3.,  0.],
       [ 0.,  3.,  0., 15.]])

In other words, the basic monomial (beyond polynomial order 1) are not orthogonal.

But if we replace the basic monomial with an explicit orthogonal polynomial constructor, we get:

hermite = chaospy.generate_expansion(3, normal)
hermite

polynomial([1.0, q0, q0**2-1.0, q0**3-3.0*q0])

hermite2 = chaospy.outer(hermite, hermite)
chaospy.E(hermite2, normal).round(15)

array([[1., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 2., 0.],
       [0., 0., 0., 6.]])

A fully diagonal matrix, which implies all the polynomials in the expansion are mutually orthogonal.

Multivariate expansions

Multivariate orthogonal expansion are (usually) created by doing a tensor product of univariate expansions together. To illustrate how this work, consider the distribution introduced in the problem formulation:

from problem_formulation import joint

joint

J(Normal(mu=1.5, sigma=0.2), Uniform(lower=0.1, upper=0.2))

Extracting the marginal density we can construct both one-dimensional expansions:

expansion0 = chaospy.generate_expansion(2, joint[0])
expansion0

polynomial([1.0, q0-1.5, q0**2-3.0*q0+2.21])

expansion1 = chaospy.generate_expansion(2, joint[1])
expansion1.round(5)

polynomial([1.0, q0-0.15, q0**2-0.3*q0+0.02167])

When constructing a multivariate expansion, it is canonical to truncate the expansion at order and graded lexicographical sorting:

chaospy.generate_expansion(2, joint).round(5)

polynomial([1.0, q1-0.15, q0-1.5, q1**2-0.3*q1+0.02167,
            q0*q1-1.5*q1-0.15*q0+0.225, q0**2-3.0*q0+2.21])

See chaospy.generate_expansion() for variations in truncations and sorting.

Wiener-Askey scheme

Polynomial chaos expansion often assume that the polynomial expansion used is of the Wiener-Askey scheme verity. The reason for this is that the expansion in the scheme correspond to orthogonality with respect to some standard probability distribution. These include:

• Hermite polynomials which are orthogonal with a normal density weight function.

• Legendre polynomials which are orthogonal with a uniform density weight function.

• Laguerre polynomials which are orthogonal with a exponential density weight function.

• Generalized Laguerre polynomials which are orthogonal with a gamma density weight function.

• Jacobi polynomials which are orthogonal with a beta density weight function.

In chaospy, these can all be constructed using chaospy.generate_expansion(). Hermite and normal distribution is showed above. The others can be created in the same way:

uniform = chaospy.Uniform(-1, 1)
legendre = chaospy.generate_expansion(3, uniform)
legendre.round(5)

polynomial([1.0, q0, q0**2-0.33333, q0**3-0.6*q0])

exponential = chaospy.Exponential()
laguerre = chaospy.generate_expansion(3, exponential)
laguerre

polynomial([1.0, q0-1.0, q0**2-4.0*q0+2.0, q0**3-9.0*q0**2+18.0*q0-6.0])

alpha = 2
gamma = chaospy.Gamma(alpha+1)
generalized_laguerre = chaospy.generate_expansion(3, gamma)
generalized_laguerre

polynomial([1.0, q0-3.0, q0**2-8.0*q0+12.0, q0**3-15.0*q0**2+60.0*q0-60.0])

alpha_, beta_ = 2, 3
beta = chaospy.Beta(alpha_+1, beta_+1, lower=-1, upper=1)
jacobi = chaospy.generate_expansion(3, beta)
jacobi.round(5)

polynomial([1.0, q0+0.14286, q0**2+0.22222*q0-0.11111,
            q0**3+0.27273*q0**2-0.27273*q0-0.0303])