Solution of a non-expensive least-square problem based on a scipy implementation

Driver:

ScipyLeastSquares

Download script:

scipy_least_squares.py

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} \]
 1import sys,os
 2import numpy as np
 3import time
 4import pandas as pd
 5import torch
 6
 7from jcmoptimizer import Server, Client, Study, Obseravtion
 8server = Server()
 9client = Client(host=server.host)
10
11
12# Definition of the search domain
13design_space = [
14    {'name': 'b1', 'type': 'continuous', 'domain': (0,10)}, 
15    {'name': 'b2', 'type': 'continuous', 'domain': (0.1,4)},
16    {'name': 'b3', 'type': 'continuous', 'domain': (-4,-0.1)},
17    {'name': 'b4', 'type': 'continuous', 'domain': (0.05,1)},
18    {'name': 'b5', 'type': 'continuous', 'domain': (0.05,1)}
19]
20constraints = [
21    {'name': 'test', 'expression': 'b2 + b3 <= 1.0'}
22]
23
24# Creation of the study object with study_id 'scipy_least_squares'
25study = client.create_study(
26    design_space=design_space,
27    constraints=constraints,
28    driver="ScipyLeastSquares",
29    study_name="Solution of a non-expensive least-square problem based on a scipy implementation",
30    study_id="scipy_least_squares"
31)
32#The vectorial model function of the MGH17 problem
33def model(x: torch.Tensor) -> torch.Tensor:
34    s = torch.arange(33)
35    return x[0] + x[1]*torch.exp(-s*x[3]) + x[2]*torch.exp(-s*x[4])
36
37#Target vector of the MGH17
38target=torch.tensor([
39    8.44E-01, 9.08E-01, 9.32E-01, 9.36E-01, 9.25E-01,
40    9.08E-01, 8.81E-01, 8.50E-01, 8.18E-01, 7.84E-01,
41    7.51E-01, 7.18E-01, 6.85E-01, 6.58E-01, 6.28E-01,
42    6.03E-01, 5.80E-01, 5.58E-01, 5.38E-01, 5.22E-01,
43    5.06E-01, 4.90E-01, 4.78E-01, 4.67E-01, 4.57E-01,
44    4.48E-01, 4.38E-01, 4.31E-01, 4.24E-01, 4.20E-01,
45    4.14E-01, 4.11E-01, 4.06E-01
46])
47
48study.configure(
49    target_vector=target.tolist(),
50    method="trf",
51    max_iter=150,
52    num_parallel=2,
53    num_initial=2,
54    jac=True
55)
56
57# Evaluation of the black-box function for specified design parameters
58def evaluate(study: Study, b1: float, b2: float, b3: float, b4: float, b5: float) -> Observation:
59
60    observation = study.new_observation()
61    #tensor of design values to reconstruct
62    x = torch.tensor([b1, b2, b3, b4, b5], requires_grad=True)
63    
64    observation.add(model(x).tolist())
65
66    #determine Jacobian matrix
67    jac = torch.autograd.functional.jacobian(
68        func=model,
69        inputs=x
70    )
71
72    for idx, param in enumerate(design_space):
73        observation.add(jac[:, idx].tolist(), derivative=param["name"])
74    return observation
75
76# Run the minimization
77study.set_evaluator(evaluate)
78study.run()
79best_sample = study.driver.best_sample
80uncertainties = study.driver.uncertainties
81print("Reconstructed parameters:")
82for param in design_space:
83    name = param['name']
84    print(f"  {name} = {best_sample[name]:.3f} +/- {uncertainties[name]:.3f}")
85