Example 2
=========

In this example, we are going to test ``cgn`` on a constrained version of the classic "Osborne problem" (taken from Wright and Holt, "Algorithms for Nonlinear Least Squares with Linear Inequality Constraints", SIAM Journal on Scientific and Statistical Computing, 1985):

$
\begin{align}
\min_{x \in \mathbb R^{11}} \quad & ||y - G(x)||_2^2 \\
\text{s. t.} \quad & C x \geq d, \\
\text{where } \quad & G_j(x) = x_1 e^{-t_j x_5} + x_2 e^{-(t_j - x_9)^2 x_6}
 + x3 e^{-(t_j - x_{10})^2 x_7} + x_4 e^{-(t_j - x_{11})^2 x_8}, \\
& t_j = \frac{j}{10}, \quad j=1,\ldots, 65, \\
& C = \left( \begin{matrix}
1 & 2 & 3 & 4 \\
1 & 1 & 0 & 0
\end{matrix} \right), \quad d = \left( \begin{matrix}
6.270063 \\ 1.741584
\end{matrix} \right)
\end{align}
$

This is a nonlinear least-squares problem with linear inequality constraints, and
its solution is known to satisfy $||y - G(x_\mathrm{min})||_2^2 = 0.0401377$.
The data vector $y$ is given by



```python
import numpy as np

y = np.array([
    1.366, 1.191, 1.112, 1.013, 0.991, 0.885, 0.831, 0.847, 0.786, 0.725,
    0.746, 0.679, 0.608, 0.655, 0.616, 0.606, 0.602, 0.626, 0.651, 0.724,
    0.649, 0.649, 0.694, 0.644, 0.624, 0.661, 0.612, 0.558, 0.533, 0.495,
    0.500, 0.423, 0.395, 0.375, 0.372, 0.391, 0.396, 0.405, 0.428, 0.429,
    0.523, 0.562, 0.607, 0.653, 0.672, 0.708, 0.633, 0.668, 0.645, 0.632,
    0.591, 0.559, 0.597, 0.625, 0.739, 0.710, 0.729, 0.720, 0.636, 0.581,
    0.428, 0.292, 0.162, 0.098, 0.054
])
```

Let us start with implementation of the function $G$ and its Jacobian.


```python
from math import exp

m = 65
n = 11

t = np.zeros(m)
for i in range(m):
    t[i] = float(i/10)

def G(x):
    z = np.ones(m)
    for i in range(m):
            z[i] = x[0]*exp(-x[4]*t[i]) + x[1]*exp(-x[5]*((t[i]-x[8])**2)) \
                + x[2]*exp(-x[6]*((t[i]-x[9])**2)) + x[3]*exp(-x[7]*((t[i]-x[10])**2))
    return z

def DG(x):
    jac = np.zeros((m, n))
    for i in range(m):
        for j in range(n):
            jac[i, 0] = exp(-x[4] * t[i])
            jac[i, 1] = exp(-x[5] * (t[i] - x[8]) ** 2)
            jac[i, 2] = exp(-x[7] * (t[i] - x[9]) ** 2)
            jac[i, 3] = exp(-x[7] * (t[i] - x[10]) ** 2)
            jac[i, 4] = - t[i] * x[0] * exp(-x[4] * t[i])
            jac[i, 5] = - (t[i] - x[8]) ** 2 * x[1] * exp(-x[5] * (t[i] - x[8]) ** 2)
            jac[i, 6] = - (t[i] - x[9]) ** 2 * x[2] * exp(-x[6] * (t[i] - x[9]) ** 2)
            jac[i, 7] = - (t[i] - x[10]) ** 2 * x[3] * exp(-x[7] * (t[i] - x[10]) ** 2)
            jac[i, 8] = 2.0 * x[5] * (t[i] - x[8]) * x[1] * exp(-x[5] * (t[i] - x[8]) ** 2)
            jac[i, 9] = 2.0 * x[6] * (t[i] - x[9]) * x[2] * exp(-x[6] * (t[i] - x[9]) ** 2)
            jac[i, 10] = 2.0 ** x[7] * (t[i] - x[10]) * x[3] * exp(-x[7] * (t[i] - x[10]) ** 2)
    return jac
```

Next, set up the matrix $C$ and the vector $d$ for the inequality constraint.



```python
C = np.zeros((2, n))
C[0, 0] = 1.0
C[0, 1] = 2.0
C[0, 3] = 3.0
C[0, 4] = 4.0
C[1, 0] = 1.0
C[1, 2] = 1.0
d = np.array([6.270063, 1.741584])
```

As initial guess, we use the canonical starting value from the cited paper:


```python
x_start = np.array([1.3, 0.65, 0.65, 0.7, 0.6, 3.0, 5.0, 7.0, 2.0, 4.5, 4.5])
```

With this, we have now all the ingredients to set up our ``cgn.Problem``.


```python
import cgn

x = cgn.Parameter(start=x_start, name="x")
incon = cgn.LinearConstraint(parameters=[x], a=C, b=d, ctype="ineq")

# Define the misfit function.
def F(x):
    return G(x) - y
# Note that DF(x) = DG(x).
problem = cgn.Problem(parameters=[x], fun=F, jac=DG, constraints=[incon])
```

To stabilize the problem, we add very mild regularization to $x$:


```python
x.beta = 1e-10
```

Now, we can solve this problem with ``cgn``:


```python
solver = cgn.CGN()
# Let the solver return some output:
solver.options.set_verbosity(lvl=2)
solution = solver.solve(problem=problem)
```

    
    
    Starting the constrained Gauss-Newton method. Cost at starting value: 3.5135292586674707
    +-----------+-------------------------+-------------------------+-------------------------+-------------------------+-------------------------+
    | Iteration | Cost                    | Constraint violation    | Stepsize (||p||)        | Steplength (h)          | Computation time [s]    |
    +-----------+-------------------------+-------------------------+-------------------------+-------------------------+-------------------------+
    +-----------+-------------------------+-------------------------+-------------------------+-------------------------+-------------------------+
    |     1     |    2.019268940326675    |           0.0           |    51.54989970703435    |          0.125          |   0.00884699821472168   |
    +-----------+-------------------------+-------------------------+-------------------------+-------------------------+-------------------------+
    |     2     |    0.6954282969801624   |           0.0           |    3.3073048651426173   |           0.5           |   0.008187055587768555  |
    +-----------+-------------------------+-------------------------+-------------------------+-------------------------+-------------------------+
    |     3     |    0.4420586847580506   |           0.0           |    8.638709316493241    |           1.0           |  0.0075435638427734375  |
    +-----------+-------------------------+-------------------------+-------------------------+-------------------------+-------------------------+
    |     4     |   0.27582424064716105   |           0.0           |    2.487560734439137    |           0.5           |   0.007577180862426758  |
    +-----------+-------------------------+-------------------------+-------------------------+-------------------------+-------------------------+
    |     5     |   0.27212108162643117   |           0.0           |    3.9421083506802677   |           0.25          |   0.00842905044555664   |
    +-----------+-------------------------+-------------------------+-------------------------+-------------------------+-------------------------+
    |     6     |   0.26265653545123335   |           0.0           |    5.363933340445492    |          0.125          |   0.008588790893554688  |
    +-----------+-------------------------+-------------------------+-------------------------+-------------------------+-------------------------+
    |     7     |    0.2462090755870182   |           0.0           |    6.0141071189406485   |          0.125          |   0.008092164993286133  |
    +-----------+-------------------------+-------------------------+-------------------------+-------------------------+-------------------------+
    |     8     |    0.2225186961624837   |           0.0           |    4.385528791791924    |          0.125          |   0.00864100456237793   |
    +-----------+-------------------------+-------------------------+-------------------------+-------------------------+-------------------------+
    |     9     |   0.19269467016533398   |           0.0           |    3.1079442522847054   |           0.25          |   0.008608818054199219  |
    +-----------+-------------------------+-------------------------+-------------------------+-------------------------+-------------------------+
    |     10    |   0.14872824443532617   |           0.0           |    2.5343909574294075   |           0.5           |   0.00801706314086914   |
    +-----------+-------------------------+-------------------------+-------------------------+-------------------------+-------------------------+
    |     11    |   0.03637188504832242   |           0.0           |    1.5181675624419728   |           1.0           |   0.007222414016723633  |
    +-----------+-------------------------+-------------------------+-------------------------+-------------------------+-------------------------+
    |     12    |   0.02043640097632607   |           0.0           |    0.5247607073097263   |           1.0           |   0.008041620254516602  |
    +-----------+-------------------------+-------------------------+-------------------------+-------------------------+-------------------------+
    |     13    |   0.020291106445352314  |           0.0           |   0.13018249760537745   |           1.0           |   0.00728607177734375   |
    +-----------+-------------------------+-------------------------+-------------------------+-------------------------+-------------------------+
    |     14    |   0.02028813024928513   |           0.0           |   0.05609490799951392   |           1.0           |   0.007730722427368164  |
    +-----------+-------------------------+-------------------------+-------------------------+-------------------------+-------------------------+
    |     15    |   0.020288130249285115  |           0.0           |   0.020203843392749696  |  9.094947017729282e-13  |   0.017717599868774414  |
    +-----------+-------------------------+-------------------------+-------------------------+-------------------------+-------------------------+
    The iteration converged successfully after 15 steps.
    Minimum of cost function: 0.020288130249285115


Let's check how our solution compares to the theoretical optimum:


```python
theoretical_minimum = 4.01377e-2
solution.cost
```




    0.020288130249285115



Wait, this is below the theoretical optimum! What's wrong? Well, recall that ``cgn`` uses the cost function
$
J(x) = \frac{1}{2} ||F(x)||_2^2 + \frac{\beta}{2} ||R(x-m)||_2^2.
$
Hence, we need to multiply this by $2$ to get the desired quantity:


```python
2 * solution.cost
```




    0.04057626049857023



This is satisfyingly close to the theoretical optimum.

01cf4486bd67292c2714efa97901e0bced58b543
ipynb
Jupyter Notebook
examples/example2.ipynb
FabianKP/cgn
9963e60c4a4bf4f3869e43d1dfbe11da74887ba5
MIT
2022-03-21T00:40:23.000Z
2022-03-21T00:40:23.000Z
examples/example2.ipynb
FabianKP/cgn
9963e60c4a4bf4f3869e43d1dfbe11da74887ba5
MIT
examples/example2.ipynb
FabianKP/cgn
9963e60c4a4bf4f3869e43d1dfbe11da74887ba5
MIT
Qwen/Qwen-72B
1. YES
2. YES