.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_surrogate_modeling/gaussian_process_regression/plot_gpr_hyperparameters_optimization.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_surrogate_modeling_gaussian_process_regression_plot_gpr_hyperparameters_optimization.py: Gaussian process fitter: configure the optimization solver ========================================================== .. GENERATED FROM PYTHON SOURCE LINES 7-44 The goal of this example is to show how to fine-tune the optimization solver used to estimate the hyperparameters of the covariance model of the Gaussian process regression. Introduction ------------ In a Gaussian process regression, there are various types of parameters which are estimated from the data. * The parameters :math:`{\bf \beta}` associated with the deterministic trend. These parameters are computed based on linear least squares. * The parameters of the covariance model. We only consider the following parameters of the covariance model: * The magnitude parameter :math:`\sigma^2` is estimated from the data. If the output dimension is equal to one, this parameter is estimated using the analytic variance estimator which maximizes the likelihood. Otherwise, if output dimension is greater than one or analytical sigma disabled, this parameter is estimated from numerical optimization. * The scale parameters :math:`{\bf \theta}\in\mathbb{R}^d` where :math:`d` is the spatial dimension of the covariance model. Often, the parameter :math:`{\bf \theta}` is a scale parameter. This step involves an optimization algorithm. All these parameters are estimated with the :class:`~openturns.experimental.GaussianProcessFitter` class. The estimation of the :math:`\vect{\theta}` parameters is the step which has the highest CPU cost. Moreover, the maximization of likelihood may be associated with difficulties e.g. many local maximums or even the non convergence of the optimization algorithm. In this case, it might be useful to fine tune the optimization algorithm so that the convergence of the optimization algorithm is, hopefully, improved. Furthermore, there are several situations in which the optimization can be initialized or completely bypassed. Suppose for example that we have already created an initial Gaussian process regression with :math:`N` points and we want to add a single new point. * It might be interesting to initialize the optimization algorithm with the optimum found for the previous Gaussian process regression: this may reduce the number of iterations required to maximize the likelihood. * We may as well completely bypass the optimization step: if the previous covariance model was correctly estimated, the update of the parameters may or may not significantly improve the estimates. This is why the goal of this example is to see how to configure the optimization of the hyperparameters of a Gaussian process regression. .. GENERATED FROM PYTHON SOURCE LINES 46-48 Definition of the model ----------------------- .. GENERATED FROM PYTHON SOURCE LINES 50-55 .. code-block:: Python import openturns as ot import openturns.experimental as otexp import openturns.viewer as otv .. GENERATED FROM PYTHON SOURCE LINES 56-57 We define the symbolic function which evaluates the output `Y` depending on the inputs `E`, `F`, `L` and `I`. .. GENERATED FROM PYTHON SOURCE LINES 59-61 .. code-block:: Python model = ot.SymbolicFunction(["E", "F", "L", "I"], ["F*L^3/(3*E*I)"]) .. GENERATED FROM PYTHON SOURCE LINES 62-63 Then we define the distribution of the input random vector. .. GENERATED FROM PYTHON SOURCE LINES 65-66 Young's modulus `E` .. GENERATED FROM PYTHON SOURCE LINES 66-79 .. code-block:: Python E = ot.Beta(0.9, 2.27, 2.5e7, 5.0e7) # in N/m^2 E.setDescription("E") # Load F F = ot.LogNormal() # in N F.setParameter(ot.LogNormalMuSigma()([30.0e3, 9e3, 15.0e3])) F.setDescription("F") # Length L L = ot.Uniform(250.0, 260.0) # in cm L.setDescription("L") # Moment of inertia I II = ot.Beta(2.5, 1.5, 310, 450) # in cm^4 II.setDescription("I") .. GENERATED FROM PYTHON SOURCE LINES 80-81 Finally, we define the dependency using a :class:`~openturns.NormalCopula`. .. GENERATED FROM PYTHON SOURCE LINES 83-89 .. code-block:: Python dim = 4 # number of inputs R = ot.CorrelationMatrix(dim) R[2, 3] = -0.2 myCopula = ot.NormalCopula(ot.NormalCopula.GetCorrelationFromSpearmanCorrelation(R)) myDistribution = ot.JointDistribution([E, F, L, II], myCopula) .. GENERATED FROM PYTHON SOURCE LINES 90-92 Create the design of experiments -------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 94-97 We consider a simple Monte-Carlo sampling as a design of experiments. This is why we generate an input sample using the :meth:`~openturns.Distribution.getSample` method of the distribution. Then we evaluate the output using the `model` function. .. GENERATED FROM PYTHON SOURCE LINES 99-103 .. code-block:: Python sampleSize_train = 10 X_train = myDistribution.getSample(sampleSize_train) Y_train = model(X_train) .. GENERATED FROM PYTHON SOURCE LINES 104-106 Create the metamodel -------------------- .. GENERATED FROM PYTHON SOURCE LINES 108-113 In order to create the Gaussian process regression, we first select a constant trend with the :class:`~openturns.LinearBasisFactory` class. Then we use a squared exponential covariance model. Finally, we use the :class:`~openturns.experimental.GaussianProcessFitter` class to calibrate a covariance model for a Gaussian process regression, taking the training sample, the covariance model and the trend basis as input arguments. The scale parameters must be positive. .. GENERATED FROM PYTHON SOURCE LINES 113-128 .. code-block:: Python dimension = myDistribution.getDimension() basis = ot.LinearBasisFactory(dimension).build() x_range = X_train.computeRange() print("x_range:") print(x_range) scale_max_factor = 4.0 # Must be > 1, tune this to match your problem scale_min_factor = 0.1 # Must be < 1, tune this to match your problem maximum_scale_bounds = scale_max_factor * x_range covarianceModel = ot.SquaredExponential([1.0] * dimension, [1.0]) algo = otexp.GaussianProcessFitter(X_train, Y_train, covarianceModel, basis) algo.run() result = algo.getResult() .. rst-class:: sphx-glr-script-out .. code-block:: none x_range: [1.31501e+07,13652.8,8.1952,63.8489] .. GENERATED FROM PYTHON SOURCE LINES 129-133 The :meth:`~openturns.experimental.GaussianProcessFitter.run` method has optimized the hyperparameters of the metamodel. We can then print the constant trend of the metamodel, which have been estimated using the least squares method. .. GENERATED FROM PYTHON SOURCE LINES 135-137 .. code-block:: Python result.getTrendCoefficients() .. raw:: html
class=Point name=Unnamed dimension=5 values=[0.000358703,-3.97972e-07,0.000413797,0.0940969,-0.026112]


.. GENERATED FROM PYTHON SOURCE LINES 138-139 We can also print the hyperparameters of the covariance model, which have been estimated by maximizing the likelihood. .. GENERATED FROM PYTHON SOURCE LINES 141-144 .. code-block:: Python basic_covariance_model = result.getCovarianceModel() print(basic_covariance_model) .. rst-class:: sphx-glr-script-out .. code-block:: none SquaredExponential(scale=[5.00377e+06,8209.1,5.80143,127.698], amplitude=[0.716676]) .. GENERATED FROM PYTHON SOURCE LINES 145-147 Get the optimizer algorithm --------------------------- .. GENERATED FROM PYTHON SOURCE LINES 149-153 The :meth:`~openturns.experimental.GaussianProcessFitter.getOptimizationAlgorithm` method returns the optimization algorithm used to optimize the :math:`\vect{\theta}` parameters of the :class:`~openturns.SquaredExponential` covariance model. .. GENERATED FROM PYTHON SOURCE LINES 155-157 .. code-block:: Python solver = algo.getOptimizationAlgorithm() .. GENERATED FROM PYTHON SOURCE LINES 158-159 Get the default optimizer. .. GENERATED FROM PYTHON SOURCE LINES 161-164 .. code-block:: Python solverImplementation = solver.getImplementation() solverImplementation.getClassName() .. rst-class:: sphx-glr-script-out .. code-block:: none 'Cobyla' .. GENERATED FROM PYTHON SOURCE LINES 165-170 The :meth:`~openturns.experimental.GaussianProcessFitter.getOptimizationBounds` method returns the bounds. The dimension of these bounds correspond to the spatial dimension of the covariance model. In the metamodeling context, this correspond to the input dimension of the model. .. GENERATED FROM PYTHON SOURCE LINES 170-173 .. code-block:: Python bounds = algo.getOptimizationBounds() bounds.getDimension() .. rst-class:: sphx-glr-script-out .. code-block:: none 4 .. GENERATED FROM PYTHON SOURCE LINES 174-178 .. code-block:: Python lbounds = bounds.getLowerBound() print("lbounds") print(lbounds) .. rst-class:: sphx-glr-script-out .. code-block:: none lbounds [13150.1,13.6528,0.0081952,0.0638489] .. GENERATED FROM PYTHON SOURCE LINES 179-183 .. code-block:: Python ubounds = bounds.getUpperBound() print("ubounds") print(ubounds) .. rst-class:: sphx-glr-script-out .. code-block:: none ubounds [2.63003e+07,27305.6,16.3904,127.698] .. GENERATED FROM PYTHON SOURCE LINES 184-185 The :meth:`~openturns.experimental.GaussianProcessFitter.getOptimizeParameters` method returns `True` if these parameters are to be optimized. .. GENERATED FROM PYTHON SOURCE LINES 185-189 .. code-block:: Python isOptimize = algo.getOptimizeParameters() print(isOptimize) .. rst-class:: sphx-glr-script-out .. code-block:: none True .. GENERATED FROM PYTHON SOURCE LINES 190-192 Configure the starting point of the optimization ------------------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 194-195 The starting point of the optimization is based on the parameters of the covariance model. .. GENERATED FROM PYTHON SOURCE LINES 195-198 .. code-block:: Python covarianceModel = ot.SquaredExponential(maximum_scale_bounds, [1.0]) algo = otexp.GaussianProcessFitter(X_train, Y_train, covarianceModel, basis) .. GENERATED FROM PYTHON SOURCE LINES 199-201 .. code-block:: Python algo.run() .. GENERATED FROM PYTHON SOURCE LINES 202-206 .. code-block:: Python result = algo.getResult() new_covariance_model = result.getCovarianceModel() print(new_covariance_model) .. rst-class:: sphx-glr-script-out .. code-block:: none SquaredExponential(scale=[5.00377e+06,8209.1,5.80143,127.698], amplitude=[0.716676]) .. GENERATED FROM PYTHON SOURCE LINES 207-208 In order to see the difference with the default optimization, we print the previous optimized covariance model. .. GENERATED FROM PYTHON SOURCE LINES 210-212 .. code-block:: Python print(basic_covariance_model) .. rst-class:: sphx-glr-script-out .. code-block:: none SquaredExponential(scale=[5.00377e+06,8209.1,5.80143,127.698], amplitude=[0.716676]) .. GENERATED FROM PYTHON SOURCE LINES 213-214 We draw the graphs that validate the meta model. .. GENERATED FROM PYTHON SOURCE LINES 216-225 .. code-block:: Python algo_gpr = otexp.GaussianProcessRegression(result) algo_gpr.run() metamodel = algo_gpr.getResult().getMetaModel() X_test = myDistribution.getSample(100) Y_test = model(X_test) validation = ot.MetaModelValidation(Y_test, metamodel(X_test)) g = validation.drawValidation() view = otv.View(g) .. image-sg:: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_hyperparameters_optimization_001.svg :alt: Metamodel validation - n = 100 :srcset: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_hyperparameters_optimization_001.svg :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 226-227 We observe that this does not change much the values of the parameters in this case. .. GENERATED FROM PYTHON SOURCE LINES 229-231 Disabling the optimization -------------------------- .. GENERATED FROM PYTHON SOURCE LINES 233-236 It is sometimes useful to completely disable the optimization of the parameters. In order to see the effect of this, we first initialize the parameters of the covariance model with the arbitrary values `[12.0, 34.0, 56.0, 78.0]`. .. GENERATED FROM PYTHON SOURCE LINES 236-239 .. code-block:: Python covarianceModel = ot.SquaredExponential([12.0, 34.0, 56.0, 78.0], [91.0]) algo = otexp.GaussianProcessFitter(X_train, Y_train, covarianceModel, basis) .. GENERATED FROM PYTHON SOURCE LINES 240-242 The :meth:`~openturns.experimental.GaussianProcessFitter.setOptimizeParameters` method can be used to disable the optimization of the parameters. .. GENERATED FROM PYTHON SOURCE LINES 242-244 .. code-block:: Python algo.setOptimizeParameters(False) .. GENERATED FROM PYTHON SOURCE LINES 245-246 Then we run the algorithm and get the result. .. GENERATED FROM PYTHON SOURCE LINES 246-249 .. code-block:: Python algo.run() result = algo.getResult() .. GENERATED FROM PYTHON SOURCE LINES 250-252 We observe that the covariance model is unchanged: the parameters have not been optimized, as required. .. GENERATED FROM PYTHON SOURCE LINES 252-255 .. code-block:: Python updated_covariance_model = result.getCovarianceModel() print(updated_covariance_model) .. rst-class:: sphx-glr-script-out .. code-block:: none SquaredExponential(scale=[12,34,56,78], amplitude=[91]) .. GENERATED FROM PYTHON SOURCE LINES 256-257 The trend, however, is still optimized, using linear least squares. .. GENERATED FROM PYTHON SOURCE LINES 257-259 .. code-block:: Python result.getTrendCoefficients() .. raw:: html
class=Point name=Unnamed dimension=5 values=[0.000319623,-3.94857e-07,0.000438513,0.0752388,-0.0166714]


.. GENERATED FROM PYTHON SOURCE LINES 260-267 Reuse the parameters from a previous optimization ------------------------------------------------- In this example, we show how to reuse the optimized parameters of a previous Gaussian process fit and configure a new one. Furthermore, we disable the optimization so that the parameters of the covariance model are not updated. This makes the process of adding a new point very fast: it improves the quality by adding a new point in the design of experiments without paying the price of the update of the covariance model. .. GENERATED FROM PYTHON SOURCE LINES 269-270 Step 1: Run a first Gaussian process fitter algorithm. .. GENERATED FROM PYTHON SOURCE LINES 270-277 .. code-block:: Python dimension = myDistribution.getDimension() basis = ot.ConstantBasisFactory(dimension).build() covarianceModel = ot.SquaredExponential(maximum_scale_bounds, [1.0]) algo = otexp.GaussianProcessFitter(X_train, Y_train, covarianceModel, basis) algo.run() result = algo.getResult() .. GENERATED FROM PYTHON SOURCE LINES 278-279 Now retrieve the optimized the covariance model. .. GENERATED FROM PYTHON SOURCE LINES 279-282 .. code-block:: Python covarianceModel = result.getCovarianceModel() print(covarianceModel) .. rst-class:: sphx-glr-script-out .. code-block:: none SquaredExponential(scale=[1.51599e+07,16372.7,16.3904,127.698], amplitude=[3.99673]) .. GENERATED FROM PYTHON SOURCE LINES 283-284 Step 2: Create a new point and add it to the previous training design. .. GENERATED FROM PYTHON SOURCE LINES 286-289 .. code-block:: Python X_new = myDistribution.getSample(20) Y_new = model(X_new) .. GENERATED FROM PYTHON SOURCE LINES 290-293 .. code-block:: Python X_train.add(X_new) X_train.getSize() .. rst-class:: sphx-glr-script-out .. code-block:: none 30 .. GENERATED FROM PYTHON SOURCE LINES 294-297 .. code-block:: Python Y_train.add(Y_new) Y_train.getSize() .. rst-class:: sphx-glr-script-out .. code-block:: none 30 .. GENERATED FROM PYTHON SOURCE LINES 298-299 Step 3: Create an updated Gaussian process fit, using the new point with the old covariance parameters. .. GENERATED FROM PYTHON SOURCE LINES 299-307 .. code-block:: Python algo = otexp.GaussianProcessFitter(X_train, Y_train, covarianceModel, basis) algo.setOptimizeParameters(False) algo.run() result = algo.getResult() notUpdatedCovarianceModel = result.getCovarianceModel() print(notUpdatedCovarianceModel) .. rst-class:: sphx-glr-script-out .. code-block:: none SquaredExponential(scale=[1.51599e+07,16372.7,16.3904,127.698], amplitude=[3.99673]) .. GENERATED FROM PYTHON SOURCE LINES 308-317 .. code-block:: Python def printCovarianceParameterChange(covarianceModel1, covarianceModel2): parameters1 = covarianceModel1.getFullParameter() parameters2 = covarianceModel2.getFullParameter() for i in range(parameters1.getDimension()): deltai = abs(parameters1[i] - parameters2[i]) print(f"Change in the parameter #{i} = {deltai}") return .. GENERATED FROM PYTHON SOURCE LINES 318-320 .. code-block:: Python printCovarianceParameterChange(covarianceModel, notUpdatedCovarianceModel) .. rst-class:: sphx-glr-script-out .. code-block:: none Change in the parameter #0 = 0.0 Change in the parameter #1 = 0.0 Change in the parameter #2 = 0.0 Change in the parameter #3 = 0.0 Change in the parameter #4 = 0.0 Change in the parameter #5 = 0.0 .. GENERATED FROM PYTHON SOURCE LINES 321-325 We see that the parameters did not change *at all*: disabling the optimization allows one to keep a constant covariance model. In a practical algorithm, we may, for example, add a block of 10 new points before updating the parameters of the covariance model. At this point, we may reuse the previous covariance model so that the optimization starts from a better point, compared to the parameters default values. This will reduce the cost of the optimization. .. GENERATED FROM PYTHON SOURCE LINES 327-329 Configure the local optimization solver --------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 331-333 The following example shows how to set the local optimization solver. We choose the `SLSQP` algorithm from :class:`~openturns.NLopt`. .. GENERATED FROM PYTHON SOURCE LINES 333-341 .. code-block:: Python problem = solver.getProblem() local_solver = ot.NLopt(problem, "LD_SLSQP") covarianceModel = ot.SquaredExponential(maximum_scale_bounds, [1.0]) algo = otexp.GaussianProcessFitter(X_train, Y_train, covarianceModel, basis) algo.setOptimizationAlgorithm(local_solver) algo.run() result = algo.getResult() .. GENERATED FROM PYTHON SOURCE LINES 342-345 .. code-block:: Python finetune_covariance_model = result.getCovarianceModel() print(finetune_covariance_model) .. rst-class:: sphx-glr-script-out .. code-block:: none SquaredExponential(scale=[2.07767e+07,48311.5,17.998,211.362], amplitude=[12.3949]) .. GENERATED FROM PYTHON SOURCE LINES 346-349 .. code-block:: Python printCovarianceParameterChange(finetune_covariance_model, basic_covariance_model) .. rst-class:: sphx-glr-script-out .. code-block:: none Change in the parameter #0 = 15772937.89251867 Change in the parameter #1 = 40102.361152787664 Change in the parameter #2 = 12.196545577041487 Change in the parameter #3 = 83.66413540622406 Change in the parameter #4 = 0.0 Change in the parameter #5 = 11.678242171996956 .. GENERATED FROM PYTHON SOURCE LINES 350-352 Configure the global optimization solver ---------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 354-361 The following example checks the robustness of the optimization of the Gaussian process fitter with respect to the optimization of the likelihood function in the covariance model parameters estimation. We use a :class:`~openturns.MultiStart` algorithm in order to avoid to be trapped by a local minimum. Furthermore, we generate the design of experiments using a :class:`~openturns.LHSExperiment`, which guarantees that the points will fill the space. .. GENERATED FROM PYTHON SOURCE LINES 363-367 .. code-block:: Python sampleSize_train = 10 X_train = myDistribution.getSample(sampleSize_train) Y_train = model(X_train) .. GENERATED FROM PYTHON SOURCE LINES 368-371 First, we create a multivariate distribution, based on independent :class:`~openturns.Uniform` marginals which have the bounds required by the covariance model. .. GENERATED FROM PYTHON SOURCE LINES 371-374 .. code-block:: Python distributions = [ot.Uniform(lbounds[i], ubounds[i]) for i in range(dim)] boundedDistribution = ot.JointDistribution(distributions) .. GENERATED FROM PYTHON SOURCE LINES 375-377 We first generate a Latin Hypercube Sampling (LHS) design made of 25 points in the sample space. This LHS is optimized so as to fill the space. .. GENERATED FROM PYTHON SOURCE LINES 379-388 .. code-block:: Python K = 25 # design size LHS = ot.LHSExperiment(boundedDistribution, K) SA_profile = ot.GeometricProfile(10.0, 0.95, 20000) LHS_optimization_algo = ot.SimulatedAnnealingLHS(LHS, ot.SpaceFillingC2(), SA_profile) LHS_optimization_algo.generate() LHS_design = LHS_optimization_algo.getResult() starting_points = LHS_design.getOptimalDesign() starting_points.getSize() .. rst-class:: sphx-glr-script-out .. code-block:: none 25 .. GENERATED FROM PYTHON SOURCE LINES 389-390 Then we create a :class:`~openturns.MultiStart` algorithm based on the LHS starting points. .. GENERATED FROM PYTHON SOURCE LINES 390-393 .. code-block:: Python solver.setMaximumIterationNumber(10000) multiStartSolver = ot.MultiStart(solver, starting_points) .. GENERATED FROM PYTHON SOURCE LINES 394-397 Finally, we configure the optimization algorithm so as to use the :class:`~openturns.MultiStart` algorithm. The bounds of the algorithm are set to match the range of the distribution used to generate the starting points. .. GENERATED FROM PYTHON SOURCE LINES 397-402 .. code-block:: Python algo = otexp.GaussianProcessFitter(X_train, Y_train, covarianceModel, basis) algo.setOptimizationAlgorithm(multiStartSolver) algo.run() result = algo.getResult() .. GENERATED FROM PYTHON SOURCE LINES 403-406 .. code-block:: Python finetune_covariance_model = result.getCovarianceModel() print(finetune_covariance_model) .. rst-class:: sphx-glr-script-out .. code-block:: none SquaredExponential(scale=[1.41989e+07,17972.2,17.0828,171.157], amplitude=[4.33969]) .. GENERATED FROM PYTHON SOURCE LINES 407-409 .. code-block:: Python printCovarianceParameterChange(finetune_covariance_model, basic_covariance_model) .. rst-class:: sphx-glr-script-out .. code-block:: none Change in the parameter #0 = 9195118.679393787 Change in the parameter #1 = 9763.078445134066 Change in the parameter #2 = 11.281394585171245 Change in the parameter #3 = 43.458991389293715 Change in the parameter #4 = 0.0 Change in the parameter #5 = 3.6230178158472963 .. GENERATED FROM PYTHON SOURCE LINES 410-411 We see that there are changes in the estimated parameters. .. GENERATED FROM PYTHON SOURCE LINES 414-415 Display figures .. GENERATED FROM PYTHON SOURCE LINES 415-416 .. code-block:: Python otv.View.ShowAll() .. _sphx_glr_download_auto_surrogate_modeling_gaussian_process_regression_plot_gpr_hyperparameters_optimization.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_gpr_hyperparameters_optimization.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_gpr_hyperparameters_optimization.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_gpr_hyperparameters_optimization.zip `