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
|
.. _polynomial_least_squares:
Least squares polynomial response surface
-----------------------------------------
Instead of replacing the model response :math:`\model(\vect{x})` with a
*local* approximation around a given set :math:`\vect{x}_0` of
input parameters as in :ref:`Taylor expansion <taylor_expansion>`, one may seek a *global* approximation of
:math:`\model(\vect{x})` over its whole domain of definition. A common
choice to this end is global polynomial approximation. For the sake of
simplicity, a *scalar* model response :math:`y=\model(\vect{x})` will
be considered from now on. Nonetheless, the following derivations hold
for a vector-valued response.
In the following, one considers global approximations of the model
response using:
- a linear function, i.e. a polynomial of degree one;
.. math::
y \approx \widehat{\model}(\vect{x}) = a_0 \, + \, \sum_{i=1}^{n_{X}} a_{i} x_i
where :math:`(a_j \, , \, j=0,\dots,n_X)` is a set of unknown
coefficients.
- a quadratic function, i.e. a polynomial of degree two.
.. math::
y \approx \widehat{\model}(\vect{x}) = a_0 \, + \, \sum_{i=1}^{n_{X}} a_{i} x_i \, + \,
\sum_{i=1}^{n_{X}} \sum_{j=1}^{n_{X}} a_{i,j} x_i x_j
The two previous equations may be recast as:
.. math::
y \approx \widehat{\model}(\vect{x}) = \sum_{j=0}^{P-1} a_j \psi_j(\vect{x})
where :math:`P` denotes the number of terms, which is equal to
:math:`n_X + 1` (resp. to :math:`1 + 2n_X + n_X (n_X - 1)/2`) when
using a linear (resp. a quadratic) approximation, and the family
:math:`(\psi_j,j=0,\dots,P-1)` gathers the constant monomial
:math:`1`, the monomials of degree one :math:`x_i` and possibly the
cross-terms :math:`x_i x_j` as well as the monomials of degree two
:math:`x_i^2`. Using the vector notation
:math:`\vect{a} = (a_{0} , \dots , a_{P-1} )^{\textsf{T}}`
and
:math:`\vect{\psi}(\vect{x}) = (\psi_{0}(\vect{x}) , \dots , \psi_{P-1}(\vect{x}) )^{\textsf{T}}`,
this can be rewritten:
.. math::
y \approx \widehat{\model}(\vect{x}) = \Tr{\vect{a}}\vect{\psi}(\vect{x})
A *global* approximation of the model response over its whole
definition domain is sought. To this end, the coefficients :math:`a_j`
may be computed using a least squares regression approach. In this
context, an experimental design, that is, a set of observations of
input parameters, is required:
.. math::
\cX = \left\{ \vect{x}^{(1)}, \dots, \vect{x}^{(N)} \right\},
as well as the corresponding model evaluations:
.. math::
\cY = \left\{ y^{(1)},\dots,y^{(N)} \right\}.
The least squares problem is to solve:
.. math::
\widehat{\vect{a}} = \argmin_{\vect{a} \in \Rset^P} \cJ(\vect{a})
where :math:`\cJ(\vect{a})` is the cost function, defined as:
.. math::
\cJ(\vect{a}) = \sum_{i=1}^N \left( y^{(i)} - \Tr{\vect{a}} \vect{\psi}\left(\vect{x}^{(i)}\right) \right)^2.
Let :math:`\vect{y} = \Tr{(y^{(1)},\dots,y^{(N)})} \in \Rset^{N}` be the
vector of output observations.
If the design matrix :math:`\mat{\Psi}` has full rank,
then the solution is given by the normal equations:
.. math::
:label: solutionLLS
\widehat{\vect{a}} = \left( \Tr{\mat{\Psi}} \mat{\Psi} \right)^{-1} \Tr{\mat{\Psi}} \vect{y}
where:
.. math::
\mat{\Psi}_{ij} = \psi_{j}\left(\vect{x}^{(i)}\right)
for :math:`i = 1, \dots, N` and :math:`j = 0, \dots, P - 1`.
A necessary condition for having a solution is that the size :math:`N`
of the experimental design is not less than the number :math:`P` of
coefficients to estimate.
The Gram matrix :math:`\Tr{\mat{\Psi}} \mat{\Psi}` can be
ill-conditionned.
Hence, the best method is not necessarily to invert the Gram matrix,
because the solution may be particularly sensitive to rounding errors.
The least-squares problem is rather solved using more robust numerical methods
such as the *singular value decomposition* (SVD) or the *QR-decomposition*.
.. topic:: API:
- See :class:`~openturns.FunctionalChaosAlgorithm`
.. topic:: Examples:
- See :doc:`/auto_meta_modeling/polynomial_chaos_metamodel/plot_functional_chaos`
.. topic:: References:
- [bjork1996]_
|
