Solution of least-square problem using Bayesian optimization

Driver:

BayesianLeastSquares

Download script: bayesian_least_squares.m

The target of the study is to showcase the solution of a non-linear least-squares problem from the NIST statistical reference datasets. As an example, the MGH17 problem is considered which consists of fitting a vectorial model \(\mathbf{f}(\mathbf{b}) \in \mathbb{R}^{33}\) with

\[f_i(\mathbf{b}) = b_1 + b_2 \exp(-i \cdot b_4) + b_3 \exp(-i \cdot b_5),\, i = 0, \dots, 32\]

to a target vector with 33 entries. The certified best-fit values are

\[ \begin{align}\begin{aligned}b_1 &= 0.3754100521 \pm 2.0723153551 \cdot 10^{-3}&\\b_2 &= 1.9358469127 \pm 0.22031669222&\\b_3 &= -1.464687136 \pm 0.22175707739&\\b_4 &= 0.1286753464 \pm 4.4861358114\cdot 10^{-3}&\\b_5 &= 0.2212269966 \pm 8.9471996575 \cdot 10^{-3}&\end{aligned}\end{align} \]

jcm_optimizer_path = '<JCM_OPTIMIZER_PATH>';
addpath(fullfile(jcm_optimizer_path, 'interface', 'matlab'));
 
server = jcmoptimizer.Server();
client = jcmoptimizer.Client(server.port); 

% Definition of the search domain
design_space = { ...
  struct('name', 'b1', 'type', 'continuous', 'domain', [0,10]), ... 
  struct('name', 'b2', 'type', 'continuous', 'domain', [0.1,4]), ...
  struct('name', 'b3', 'type', 'continuous', 'domain', [-4,-0.1]), ...
  struct('name', 'b4', 'type', 'continuous', 'domain', [0.05,1]), ...
  struct('name', 'b5', 'type', 'continuous', 'domain', [0.05,1]) ...
};
constraints = { ...
    struct('name', 'test', 'expression', 'b2 + b3 <= 1.0') ...
};

 % Creation of the study object with study_id 'bayesian_least_squares'
study = client.create_study( ...
    'design_space', design_space, ...
    'constraints', constraints,...
    'driver','BayesianLeastSquares',...
    'name','Solution of least-square problem using Bayesian optimization',...
    'study_id', 'bayesian_least_squares');
%The vectorial model function of the MGH17 problem
function val = model(x)
    s = 0:32;
    val = x(1) + x(2)*exp(-s.*x(4)) + x(3)*exp(-s.*x(5));
end

%Target vector of the MGH17
target=[8.44E-01, 9.08E-01, 9.32E-01, 9.36E-01, 9.25E-01, ...
 9.08E-01, 8.81E-01, 8.50E-01, 8.18E-01, 7.84E-01, ...
 7.51E-01, 7.18E-01, 6.85E-01, 6.58E-01, 6.28E-01, ...
 6.03E-01, 5.80E-01, 5.58E-01, 5.38E-01, 5.22E-01, ...
 5.06E-01, 4.90E-01, 4.78E-01, 4.67E-01, 4.57E-01, ...
 4.48E-01, 4.38E-01, 4.31E-01, 4.24E-01, 4.20E-01, ...
 4.14E-01, 4.11E-01, 4.06E-01];

study.configure( ...
    'target_vector', target, ...
    'max_iter', 120 ...
);
% Evaluation of the black-box function for specified design parameters
function observation = evaluate(study, sample)

    observation = study.new_observation();
    %array of design values
    x = [sample.b1, sample.b2, sample.b3, sample.b4, sample.b5];    
    observation.add(model(x));
    
end  

% Run the minimization
study.set_evaluator(@evaluate);
study.run(); 

best_sample = study.driver.best_sample;
min_chisq = study.driver.min_objective;
uncertainties = study.driver.uncertainties;
fprintf("Reconstructed parameters with chi-squared value %e\n", min_chisq);
fns = fieldnames(best_sample);
for i = 1:length(fns)
    fprintf("  %s = %f +/- %f\n", fns{i}, ...
            best_sample.(fns{i}), uncertainties.(fns{i})); 
end
 
% Before running a Markov-chain Monte-Carlo (MCMC) sampling we converge the surrogate
% models by sampling around the minimum. To make the study more explorative, the
% scaling parameter is increased and the effective degrees of freedom is set to one.
study.configure( ...
    'scaling', 10.0, ...
    'effective_DOF', 1.0,  ...
    'min_uncertainty', min_chisq*1e-8, ...
    'max_iter', 150, ...
    'min_val', 0.0 ...
);
while(not(study.is_done))
   sug = study.get_suggestion();
   obs = evaluate(study, sug.sample);
   study.add_observation(obs, sug.id);
end

% Run the MCMC sampling with 32 walkers
num_walkers = 32; max_iter = 10000;
mcmc_result = study.driver.run_mcmc( ...
    'rel_error', 0.01, ...
    'num_walkers', num_walkers, ...
    'max_iter', max_iter ...
);
% This Matlab code does not contain a detailed analysis of the MCMC sampling by corner plots.
% Please, see the corresponding Python tutorial for a code example.

Markov-Chain Monte-Carlo (MCMC) sampling of the probability density of the parameters \(b_1,\dots,b_5\) based on the analytic model function.

Markov-Chain Monte-Carlo (MCMC) sampling of the probability density of the parameters \(b_1,\dots,b_5\) based on the trained surrogate of the study. A comparison between the analytic and the surrogate model function shows a good quantitative agreement.