1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
|
.. _orthogonal_polynomials:
Orthogonal polynomials
----------------------
This page provides mathematical details on sequences of
orthogonal polynomials. Some of these sequences will be used to
construct the basis of the so-called *polynomial chaos expansion*.
Mathematical framework
~~~~~~~~~~~~~~~~~~~~~~
Orthogonal polynomials are associated to an inner product, defined
as follows.
Given an *interval of orthogonality* :math:`[\alpha,\beta]`
(:math:`\alpha \in \Rset \cup \{-\infty\}`,
:math:`\beta \in \Rset \cup \{\infty\}`, :math:`\alpha < \beta`) and a
weight function :math:`w(x)> 0`, the polynomials :math:`P`
and :math:`Q` are orthogonal if:
.. math::
\scalarproduct{P}{Q} = \int_{\alpha}^{\beta}P(x)Q(x)~w(x) dx = 0
Therefore, a sequence of orthogonal polynomials :math:`(P_n)_{n\geq 0}`
(:math:`P_n`: polynomial of degree :math:`n`) verifies:
.. math::
\scalarproduct{P_m}{P_n} = 0 \text{~~for every~~} m \neq n
The chosen inner product induces a norm on polynomials in the usual
way:
.. math::
\parallel P_n\parallel = \scalarproduct{P_n}{P_n}^{1/2}
In the following, we consider weight functions :math:`w(x)`
corresponding to *probability density functions*, which satisfy:
.. math::
\int_{\alpha}^{\beta} \; w(x) \; dx \, \, = \,\, 1
Moreover, we consider families of *orthonormal* polynomials
:math:`(P_n)_{n\geq 0}`, that is polynomials with a unit norm:
.. math::
\|P_n\| \, \, = \, \, 1
Any sequence of orthogonal polynomials has a recurrence formula
relating any three consecutive polynomials as follows:
.. math::
P_{n+1}\ =\ (a_nx+b_n)\ P_n\ +\ c_n\ P_{n-1}
Orthogonormal polynomials with respect to usual probability distributions
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Below is a table showing an example of particular (normalized)
orthogonal polynomials associated with *continuous* weight functions.
Note that the orthonormal polynomials are
orthonormal with respect to the standard representative distribution
of the given distribution.
+-----------------+------------------------------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
| Ortho. poly. | :math:`P_n(x)` | Weight :math:`w(x)^{\strut}` | Recurrence coefficients :math:`(a_n,b_n,c_n)` |
+=================+==========================================================================================+=====================================================================================================================+===================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================+
| Hermite | :math:`{He}_n(x)^{\strut}` | :math:`\displaystyle \frac{1}{\sqrt{2 \pi}} e^{-\frac{x^2}{2}}` | :math:`\begin{array}{ccc} a_n & = & \frac{1}{\sqrt{n+1}} \\ b_n & = & 0 \\ c_n & = & - \sqrt{\frac{n}{n+1}} \end{array}` |
+-----------------+------------------------------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
| Legendre | :math:`\begin{array}{c} {Le}_n(x) \\ \\ \alpha>-1 \\ \end{array}` | :math:`\displaystyle \frac{1}{2}^{\strut} \times \mathbb{I}_{[-1,1]}(x)` | :math:`\begin{array}{ccc} a_n & = & \frac{\sqrt{(2n+1)(2n+3)}}{n+1} \\ b_n & = & 0 \\ c_n & = & -\frac{ n \sqrt{2n+3} }{ (n+1)\sqrt{2n-1} } \end{array}` |
+-----------------+------------------------------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
| Laguerre | :math:`L_n^{(\alpha)}(x)` | :math:`\displaystyle \frac{x^{k-1}}{\Gamma(k)}~e^{-x} \mathbb{I}_{[0,+\infty[}(x)` | :math:`\begin{array}{ccc} \omega_{n} & = & \left((n+1)(n+k+1) \right)^{-1/2} \\ a_n & = & \omega_{n} \\ b_n & = & -(2n+k+1)~\omega_{n} \\ c_n & = & -\sqrt{(n+k)n}~\omega_{n} \end{array}` |
+-----------------+------------------------------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
| Jacobi | :math:`\begin{array}{c} J^{(\alpha,\beta)}_n(x) \\ \\ \\ \alpha,\beta>-1 \\ \end{array}` | :math:`\frac{(1-x)^{\alpha}(1+x)^{\beta}}{2^{\alpha + \beta + 1} B(\beta + 1, \alpha + 1)} \mathbb{I}_{[-1,1]}(x)` | :math:`\begin{array}{ccc} K_{1,n} & = & \frac{2n+\alpha + \beta + 3}{(n+1)(n+\alpha+1)(n+\beta+1)(n+\alpha+\beta+1)} \\ \\ K_{2,n} & = & \frac{1}{2} \sqrt{(2n + \alpha + \beta + 1) K_{1,n}} \\ \\a_n & = & K_{2,n}(2n+\alpha + \beta + 2) \\ \\ b_n & = & K_{2,n}\frac{(\alpha - \beta)(\alpha + \beta)}{2n+\alpha+\beta} \\ \\ c_n & = & - \frac{2n+\alpha+\beta + 2}{2n+\alpha+\beta} \Big[(n+\alpha)(n+\beta) \\ & & \times (n+\alpha+\beta)n\frac{K_{1,n}}{2n+\alpha+\beta-1}\Big]^{1/2} \end{array}` |
+-----------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
Furthermore, two families of orthonormal polynomials with respect to
*discrete* probability distributions are available, namely
the so-called Charlier and Krawtchouk polynomials:
+----------------------------------+-------------------------------------------------------------------------------------------+----------------------------------------------------------------------------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
| Ortho. poly. | :math:`P_n(x)` | Probability mass function | Recurrence coefficients :math:`(a_n,b_n,c_n)` |
+==================================+===========================================================================================+================================================================================================================+=========================================================================================================================================================================================================================+
| Charlier | :math:`\begin{array}{c} Ch^{(\lambda)}_n(x) \\ \\ \lambda>0 \\ \end{array}` | :math:`\begin{array}{c} \displaystyle{\frac{\lambda^k}{k!}~e^{-\lambda}} \\ \\ k=0,1,2,\dots \\ \end{array}` | :math:`\begin{array}{ccc} a_n & = & - \frac{1}{\sqrt{\lambda (n+1)}} \\ \\ b_n & = & \frac{n+\lambda}{\sqrt{\lambda (n+1)}} \\ \\ c_n & = & - \sqrt{1 - \frac{1}{n+1}} \end{array}` |
+----------------------------------+-------------------------------------------------------------------------------------------+----------------------------------------------------------------------------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
| Krawtchouk\ :math:`^{\dagger}` | :math:`\begin{array}{c} Kr^{(m,p)}_n(x) \\ \\ m \in \Nset~,~p \in [0,1] \\ \end{array}` | :math:`\begin{array}{c} \displaystyle{\binom{m}{k}p^k (1-p)^{m-k}} \\ \\ k=0,1,2,\dots \\ \end{array}` | :math:`\begin{array}{ccc} a_n & = & - \frac{1}{\sqrt{(n+1)(m-n)p(1-p)}} \\ \\ b_n & = & \frac{p(m-n)+n(1-p)}{\sqrt{(n+1)(m-n)p(1-p)}} \\ \\ c_n & = & - \sqrt{(1 - \frac{1}{n+1})(1+\frac{1}{m-n})} \end{array}` |
+----------------------------------+-------------------------------------------------------------------------------------------+----------------------------------------------------------------------------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
Notice that the Krawtchouk polynomials are only defined up to a
degree :math:`n` equal to :math:`m-1`. Indeed, for :math:`n=m`, some
factors of the denominators of the recurrence coefficients would be
equal to zero.
To sum up, the distribution type are reported in
the table below together with the associated families of orthonormal
polynomials.
+----------------------------------------+-------------------------+----------------------------------+---------------------------------------+
| Distribution | Support | Polynomial family | In the library |
+========================================+=========================+==================================+=======================================+
| Normal :math:`\cN(0,1)` | :math:`\Rset` | Hermite | :class:`~openturns.HermiteFactory` |
+----------------------------------------+-------------------------+----------------------------------+---------------------------------------+
| Uniform :math:`\cU(-1,1)` | :math:`[-1,1]` | Legendre | :class:`~openturns.LegendreFactory` |
+----------------------------------------+-------------------------+----------------------------------+---------------------------------------+
| Gamma :math:`\Gamma(k,1,0)` | :math:`(0,+\infty)` | Laguerre | :class:`~openturns.LaguerreFactory` |
+----------------------------------------+-------------------------+----------------------------------+---------------------------------------+
| Beta :math:`B(\alpha,\beta,-1,1)` | :math:`(-1,1)` | Jacobi | :class:`~openturns.JacobiFactory` |
+----------------------------------------+-------------------------+----------------------------------+---------------------------------------+
| Poisson :math:`\cP(\lambda)` | :math:`\Nset` | Charlier | :class:`~openturns.CharlierFactory` |
+----------------------------------------+-------------------------+----------------------------------+---------------------------------------+
| Binomial :math:`\cB(m,p)` | :math:`\{0,\dots,m\}` | Krawtchouk\ :math:`^{\dagger}` | :class:`~openturns.KrawtchoukFactory` |
+----------------------------------------+-------------------------+----------------------------------+---------------------------------------+
| Negative Binomial :math:`\cN \cB(m,p)` | :math:`\Nset` | Meixner | :class:`~openturns.MeixnerFactory` |
+----------------------------------------+-------------------------+----------------------------------+---------------------------------------+
Orthogonal polynomials with respect to arbitrary probability distributions
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
It is also possible to generate a family of orthonormal polynomials
with respect to an arbitrary probability distribution :math:`w(x)`.
The *Gram-Schmidt* algorithm can be used to this end. Note
that this algorithm gives a constructive proof of the existence of
orthonormal bases.
However it is known to be numerically unstable, so alternative
procedures are often used in practice. The available orthonormalization
algorithm is the *Stieltjes* algorithm.
.. topic:: API:
- See the available :ref:`orthogonal basis <orthogonal_basis>`.
.. topic:: Examples:
- See :doc:`/auto_meta_modeling/polynomial_chaos_metamodel/plot_functional_chaos`
.. topic:: References:
- [gautschi2004]_
- [chihara1978]_
- [sullivan2015]_ chapter 8 page 133
|
