Least squares polynomial response surface

Instead of replacing the model response h(\underline{x}) for a local approximation around a given set \underline{x}_0 of input parameters, one may seek a global approximation of h(\underline{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 y=h(\underline{x}) will be considered from now on. Nonetheless, the following derivations hold for a vector-valued response.
In the sequel, one considers global approximations of the model response using:
  • a linear function, i.e. a polynomial of degree one;

  • a quadratic function, i.e. a polynomial of degree two.

y \, \, \approx \, \, \widehat{h}(\underline{x}) \, \, = \, \, a_0 \, + \,  \sum_{i=1}^{n_{X}} \; a_{i} \; x_i

where (a_j  \, , \, j=0,\dots,n_X) is a set of unknown coefficients.

\begin{aligned}
    \underline{y} \, \, \approx \, \, \widehat{h}(\underline{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
  \end{aligned}

The two previous equations may be recast using a unique formalism as follows:

\underline{y} \, \, \approx \, \, \widehat{h}(\underline{x}) \, \, = \, \, \sum_{j=0}^{P-1} \; a_j \; \psi_j(\underline{x})

where P denotes the number of terms, which is equal to (n_X + 1) (resp. to (1 + 2n_X + n_X (n_X - 1)/2)) when using a linear (resp. a quadratic) approximation, and the family (\psi_j,j=0,\dots,P-1) gathers the constant monomial 1, the monomials of degree one x_i and possibly the cross-terms x_i x_j as well as the monomials of degree two x_i^2. Using the vector notation \underline{a} \, \, = \, \, (a_{0} , \dots , a_{P-1} )^{\textsf{T}} and \underline{\psi}(\underline{x}) \, \, = \, \, (\psi_{0}(\underline{x}) , \dots , \psi_{P-1}(\underline{x}) )^{\textsf{T}}, this rewrites:

\underline{y} \, \, \approx \, \, \widehat{h}(\underline{x}) \, \, = \, \, \underline{a}^{\textsf{T}} \; \underline{\psi}(\underline{x})

A global approximation of the model response over its whole definition domain is sought. To this end, the coefficients a_j may be computed using a least squares regression approach. In this context, an experimental design \underline{\cX} =(x^{(1)},\dots,x^{(N)}), i.e. a set of realizations of input parameters is required, as well as the corresponding model evaluations \underline{\cY} =(y^{(1)},\dots,y^{(N)}).
The following minimization problem has to be solved:

\mbox{Find} \quad \widehat{\underline{a}} \quad \mbox{that minimizes} \quad \cJ(\underline{a}) \, \, = \, \, \sum_{i=1}^N \; \left( y^{(i)} \; - \; \underline{a}^{\textsf{T}}  \underline{\psi}(\underline{x}^{(i)}) \right)^2

The solution is given by:

\widehat{\underline{a}} \, \, = \, \, \left( \underline{\underline{\Psi}}^{\textsf{T}} \underline{\underline{\Psi}}  \right)^{-1} \; \underline{\underline{\Psi}}^{\textsf{T}}  \; \underline{\cY}

where:

\underline{\underline{\Psi}} \, \, = \, \, (\psi_{j}(\underline{x}^{(i)}) \; , \; i=1,\dots,N \; , \; j = 0,\dots,P-1)

It is clear that the above equation is only valid for a full rank information matrix. A necessary condition is that the size N of the experimental design is not less than the number P of PC coefficients to estimate. In practice, it is not recommended to directly invert \underline{\underline{\Psi}}^{\textsf{T}} \underline{\underline{\Psi}} since the solution may be particularly sensitive to an ill-conditioning of the matrix. The least-square problem is rather solved using more robust numerical methods such as singular value decomposition (SVD) or QR-decomposition.

References:

    1. Bjorck, 1996, “Numerical methods for least squares problems”, SIAM Press, Philadelphia, PA.