Blurb:: Uncertainty quantification with stochastic collocation

Description::
Stochastic collocation is a general framework for approximate
representation of random response functions in terms of
finite-dimensional interpolation bases.  

The stochastic collocation (SC) method is very similar to 
\ref method-polynomial_chaos, with the key difference that the orthogonal
polynomial basis functions are replaced with interpolation polynomial
bases. The interpolation polynomials may be either local or global
and either value-based or gradient-enhanced. In the local case,
valued-based are piecewise linear splines and gradient-enhanced are
piecewise cubic splines, and in the global case, valued-based are
Lagrange interpolants and gradient-enhanced are Hermite interpolants. 
A value-based expansion takes the form

\f[R = \sum_{i=1}^{N_p} r_i L_i(\xi) \f]

where \f$N_p\f$ is the total number of collocation points, \f$r_i\f$
is a response value at the \f$i^{th}\f$ collocation point, \f$L_i\f$
is the \f$i^{th}\f$ multidimensional interpolation polynomial, and
\f$\xi\f$ is a vector of standardized random variables.

Thus, in PCE, one forms coefficients for known orthogonal polynomial
basis functions, whereas SC forms multidimensional interpolation
functions for known coefficients.

<!--
The following provides details on the various stochastic collocation
method options in Dakota.

The groups and optional keywords relating to method specification are: 
\li Group 1 
\li Group 2
\li Group 3
\li Group 4
\li \c dimension_preference
\li \c use_derivatives
\li \c fixed_seed

In addition, this method treats variables that are not
aleatoric-uncertain different, despite the \ref variables-active keyword.

Group 5 and the remainder of the optional keywords relate to the output 
of the method.
-->

<b> Basis polynomial family (Group 2) </b> 

In addition to the \ref method-stoch_collocation-askey and \ref
method-stoch_collocation-wiener basis types also supported by \ref
method-polynomial_chaos, SC supports the option of \c piecewise local 
basis functions. These are piecewise linear splines, or in the case of
gradient-enhanced interpolation via the \c use_derivatives
specification, piecewise cubic Hermite splines. Both of these basis
options provide local support only over the range from the
interpolated point to its nearest 1D neighbors (within a tensor grid
or within each of the tensor grids underlying a sparse grid), which
exchanges the fast convergence of global bases for smooth functions
for robustness in the representation of nonsmooth response functions
(that can induce Gibbs oscillations when using high-order global basis
functions). When local basis functions are used, the usage of
nonequidistant collocation points (e.g., the Gauss point selections
described above) is not well motivated, so equidistant Newton-Cotes
points are employed in this case, and all random variable types are
transformed to standard uniform probability space. The
global gradient-enhanced interpolants (Hermite interpolation
polynomials) are also restricted to uniform or transformed uniform
random variables (due to the need to compute collocation weights by
integration of the basis polynomials) and share the variable support
shown in \ref topic-variable_support for Piecewise SE. Due to numerical 
instability in these high-order basis polynomials, they are deactivated 
by default but can be activated by developers using a compile-time switch.

<b> Interpolation grid type (Group 3) </b> 

To form the multidimensional interpolants \f$L_i\f$ of the expansion,
two options are provided.

<ol>
<li> interpolation on a tensor-product of Gaussian quadrature points
 (specified with \c quadrature_order and, optionally, \c
 dimension_preference for anisotropic tensor grids). As for PCE,
 non-nested Gauss rules are employed by default, although the
 presence of \c p_refinement or \c h_refinement will result in
 default usage of nested rules for normal or uniform variables
 after any variable transformations have been applied (both
 defaults can be overridden using explicit \c nested or \c
 non_nested specifications).
<li> interpolation on a Smolyak sparse grid (specified with \c
 sparse_grid_level and, optionally, \c dimension_preference for
 anisotropic sparse grids) defined from Gaussian rules. As for
 sparse PCE, nested rules are employed unless overridden with the
 \c non_nested option, and the growth rules are restricted unless
 overridden by the \c unrestricted keyword.
</ol>

Another distinguishing characteristic of stochastic collocation
relative to \ref method-polynomial_chaos is the ability to reformulate the
interpolation problem from a \c nodal interpolation approach into a \c
hierarchical formulation in which each new level of interpolation
defines a set of incremental refinements (known as hierarchical
surpluses) layered on top of the interpolants from previous levels.
This formulation lends itself naturally to uniform or adaptive
refinement strategies, since the hierarchical surpluses can be
interpreted as error estimates for the interpolant. Either global or
local/piecewise interpolants in either value-based or
gradient-enhanced approaches can be formulated using \c hierarchical
interpolation. The primary restriction for the hierarchical case is
that it currently requires a sparse grid approach using nested
quadrature rules (Genz-Keister, Gauss-Patterson, or Newton-Cotes for
standard normals and standard uniforms in a transformed space: Askey,
Wiener, or Piecewise settings may be required), although this
restriction can be relaxed in the future. A selection of \c
hierarchical interpolation will provide greater precision in the
increments to mean, standard deviation, covariance, and
reliability-based level mappings induced by a grid change within
uniform or goal-oriented adaptive refinement approaches (see following
section).

It is important to note that, while \c quadrature_order and \c
sparse_grid_level are array inputs, only one scalar from these arrays
is active at a time for a particular expansion estimation.  These
scalars can be augmented with a \c dimension_preference to support
anisotropy across the random dimension set.  The array inputs are
present to support advanced use cases such as multifidelity UQ, where
multiple grid resolutions can be employed.

<b> Automated refinement type (Group 1) </b> 

Automated expansion refinement can be selected as either \c
p_refinement or \c h_refinement, and either refinement specification
can be either \c uniform or \c dimension_adaptive. The \c
dimension_adaptive case can be further specified as either \c sobol or
\c generalized (\c decay not supported). Each of these automated
refinement approaches makes use of the \c max_iterations and \c
convergence_tolerance iteration controls.
The \c h_refinement specification involves use of the same piecewise
interpolants (linear or cubic Hermite splines) described above for the
\c piecewise specification option (it is not necessary to redundantly
specify \c piecewise in the case of \c h_refinement). In future
releases, the \c hierarchical interpolation approach will enable local
refinement in addition to the current \c uniform and \c
dimension_adaptive options.

<b> Covariance type (Group 5) </b> 

These two keywords are used to specify how this method computes, stores,
and outputs the covariance of the responses.  In particular, the diagonal
covariance option is provided for reducing post-processing overhead and 
output volume in high dimensional applications.

<b> Active Variables </b>

The default behavior is to form expansions over aleatory 
uncertain continuous variables. To form expansions 
over a broader set of variables, one needs to specify 
\c active followed by \c state, \c epistemic, \c design, or \c all 
in the variables specification block. 

For continuous design, continuous state, and continuous epistemic
uncertain variables included in the expansion, interpolation points
for these dimensions are based on Gauss-Legendre rules if non-nested,
Gauss-Patterson rules if nested, and Newton-Cotes points in the case
of piecewise bases. Again, when probability integrals are evaluated,
only the aleatory random variable domain is integrated, leaving behind
a polynomial relationship between the statistics and the remaining
design/state/epistemic variables.

<b> Optional Keywords regarding method outputs </b>

Each of these sampling specifications refer to sampling on the SC
approximation for the purposes of generating approximate statistics.
\li \c sample_type
\li \c samples
\li \c seed
\li \c fixed_seed
\li \c rng
\li \c probability_refinement 
\li \c distribution
\li \c reliability_levels
\li \c response_levels
\li \c probability_levels
\li \c gen_reliability_levels

Since SC approximations are formed on structured grids, there should
be no ambiguity with simulation sampling for generating the SC expansion.

When using the \c probability_refinement control, the number of
refinement samples is not under the user's control (these evaluations
are approximation-based, so management of this expense is less
critical). This option allows for refinement of probability and
generalized reliability results using importance sampling.

<b> Multi-fidelity UQ </b>

When using multifidelity UQ, the high fidelity expansion generated
from combining the low fidelity and discrepancy expansions retains the
polynomial form of the low fidelity expansion (only the coefficients
are updated).  Refer to \ref method-polynomial_chaos for information
on the multifidelity interpretation of array inputs for \c
quadrature_order and \c sparse_grid_level.

<b> Usage Tips </b>

If \e n is small, then tensor-product Gaussian quadrature is again the
preferred choice. For larger \e n, tensor-product quadrature quickly
becomes too expensive and the sparse grid approach is preferred. For
self-consistency in growth rates, nested rules employ restricted
exponential growth (with the exception of the \c dimension_adaptive \c
p_refinement \c generalized case) for consistency with the linear
growth used for non-nested Gauss rules (integrand precision
\f$i=4l+1\f$ for sparse grid level \e l and \f$i=2m-1\f$ for tensor
grid order \e m).

<b> Additional Resources </b>

%Dakota provides access to SC methods through the NonDStochCollocation
class. Refer to the Uncertainty Quantification Capabilities chapter of
the Users Manual \cite UsersMan and the Stochastic Expansion Methods
chapter of the Theory Manual \cite TheoMan for additional information
on the SC algorithm.

<b> Expected HDF5 Output </b>

If Dakota was built with HDF5 support and run with the 
\ref environment-results_output-hdf5 keyword, this method
writes the following results to HDF5:

- \ref hdf5_results-se_moments
- \ref hdf5_results-pdf
- \ref hdf5_results-level_mappings
- \ref hdf5_results-vbd


Topics::

Examples::
\verbatim
method,
	stoch_collocation
	  sparse_grid_level = 2	
	  samples = 10000 seed = 12347 rng rnum2	
	  response_levels = .1 1. 50. 100. 500. 1000.	
	  variance_based_decomp
\endverbatim

Theory::
As mentioned above, a value-based expansion takes the form

\f[R = \sum_{i=1}^{N_p} r_i L_i(\xi) \f]

The \f$i^{th}\f$ interpolation polynomial assumes the value of 1 at
the \f$i^{th}\f$ collocation point and 0 at all other collocation
points, involving either a global Lagrange polynomial basis or local
piecewise splines. It is easy to see that the approximation reproduces
the response values at the collocation points and interpolates between
these values at other points. A gradient-enhanced expansion (selected
via the \c use_derivatives keyword) involves both type 1 and type 2
basis functions as follows:

\f[R = \sum_{i=1}^{N_p} [ r_i H^{(1)}_i(\xi)
 + \sum_{j=1}^n \frac{dr_i}{d\xi_j} H^{(2)}_{ij}(\xi) ] \f]

where the \f$i^{th}\f$ type 1 interpolant produces 1 for the value at
the \f$i^{th}\f$ collocation point, 0 for values at all other
collocation points, and 0 for derivatives (when differentiated) at all
collocation points, and the \f$ij^{th}\f$ type 2 interpolant produces
0 for values at all collocation points, 1 for the \f$j^{th}\f$
derivative component at the \f$i^{th}\f$ collocation point, and 0 for
the \f$j^{th}\f$ derivative component at all other collocation points.
Again, this expansion reproduces the response values at each of the
collocation points, and when differentiated, also reproduces each
component of the gradient at each of the collocation points. Since
this technique includes the derivative interpolation explicitly, it
eliminates issues with matrix ill-conditioning that can occur in the
gradient-enhanced PCE approach based on regression. However, the
calculation of high-order global polynomials with the desired
interpolation properties can be similarly numerically challenging such
that the use of local cubic splines is recommended due to numerical
stability.

<!-- Rhe orthogonal polynomials used in defining
the Gauss points that make up the interpolation grid are governed by
the Wiener, Askey, or Extended options. The Wiener option uses
interpolation points from Gauss-Hermite (non-nested) or Genz-Keister
(nested) integration rules for all random variables and employs the
same nonlinear variable transformation as the local and global
reliability methods (and therefore has the same variable support).
The Askey option, however, employs interpolation points from
Gauss-Hermite (Genz-Keister if nested), Gauss-Legendre
(Gauss-Patterson if nested), Gauss-Laguerre, Gauss-Jacobi, and
generalized Gauss-Laguerre quadrature. The Extended option avoids the
use of any nonlinear variable transformations by augmenting the Askey
approach with Gauss points from numerically-generated orthogonal
polynomials for non-Askey probability density functions. As for PCE,
the Wiener/Askey/Extended selection defaults to Extended, can be
overridden by the user using the keywords \c askey or \c wiener, and
automatically falls back from Extended/Askey to Wiener on a per
variable basis as needed to support prescribed correlations.-->

Faq::
See_Also::	method-adaptive_sampling, method-gpais, method-local_reliability, method-global_reliability, method-sampling, method-importance_sampling, method-polynomial_chaos
