Example 1
=========

Let us test ``cgn`` on a simple unconstrained linear least-squares problem. As an example, we use a regularized version of problem 32 from the article

More, J., Garbox, B. and Hillstrom, E. "Testing Unconstrained Optimization Software", 1981.

The problem is as follows:

$
\begin{align}
\min_{x \in \mathbb R^n} & \frac{1}{2}||F(x)||_2^2 + \frac{\beta}{2}||x||_2^2, \\
\text{where} \quad & F(x)_j = x_j - \frac{2}{m} \sum_{i=1}^n x_i - 1, \quad \text{for } 1 \leq j \leq n, \\
& F(x)_j = - \frac{2}{m} \sum_{i=1}^n x_i - 1, \quad \text{for } n < j <= m, \\
& \beta = \frac{1}{100},
\end{align}
$

with $m >= n$. Let us choose $n=200$ and $m=400$.

First, we implement the affine misfit function $F=F(x)$.


```python
import numpy as np

m = 400
n = 200

def F(x):
    z = np.zeros(m)
    for i in range(n):
        z[i] = x[i] - 2.0 * sum(x) / m - 1.
    for i in range(n,m):
        z[i] = - 2.0 * sum(x) / m - 1.
    return z
```

We also have to implement its Jacobian:



```python
A = np.zeros((m,n))
# upper half of the Jacobian matrix
for i in range(n):
    for j in range(n):
        if i == j:
            A[i,j] = 1.0 - 2.0 / m
        else:
            A[i,j] = -2.0 / m
# lower half of the Jacobian matrix
for i in  range(n,m):
    for j in range(n):
        A[i,j] = -2.0 / m

def DF(x):
    return A
```

Let us now set up the ``cgn.Parameter`` object. To create a ``cgn.Parameter`` object, the user has to provide a starting guess at initialization. This is necessary so that ``cgn`` can verify the consistency of the user input (e.g. that the dimensions match) before it starts solving the optimization problem. By doing that, it makes it easier to diagnose runtime errors that occur later. For our parameter $x$, we simply use the 1-vector as initial guess.


```python
import cgn

x = cgn.Parameter(start=np.ones(n), name="x")
```

Note that the starting guess can be adapted later via the ``x.start``-property.
Next, we set up the regularization term by setting $\beta = 1$. The regularization operator defaults to the identity, and the regularizing guess default to the zero vector.


```python
x.beta = 0.001
```

Next, we can initialize the ``cgn.Problem`` object:


```python
problem = cgn.Problem(parameters=[x], fun=F, jac=DF)
```

We initialize the solver...


```python
solver = cgn.CGN()
```

and solve the problem:


```python
x_start = np.ones(n)
solver.options.set_verbosity(2)
solution = solver.solve(problem=problem)
```

    
    
    Starting the constrained Gauss-Newton method. Cost at starting value: 500.1
    +-----------+-------------------------+-------------------------+-------------------------+-------------------------+-------------------------+
    | Iteration | Cost                    | Constraint violation    | Stepsize (||p||)        | Steplength (h)          | Computation time [s]    |
    +-----------+-------------------------+-------------------------+-------------------------+-------------------------+-------------------------+
    +-----------+-------------------------+-------------------------+-------------------------+-------------------------+-------------------------+
    |     1     |    100.09990009990057   |           0.0           |    28.270143239845737   |           1.0           |   0.012190103530883789  |
    +-----------+-------------------------+-------------------------+-------------------------+-------------------------+-------------------------+
    |     2     |    100.0999000999005    |           0.0           |   1.55683666436841e-13  |           0.25          |   0.03589987754821777   |
    +-----------+-------------------------+-------------------------+-------------------------+-------------------------+-------------------------+
    The iteration converged successfully after 2 steps.
    Minimum of cost function: 100.0999000999005


Let us view the solution:


```python
x_min = solution.minimizer("x")
x_min
```




    array([-0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001, -0.999001, -0.999001, -0.999001, -0.999001,
           -0.999001, -0.999001])



Let us compare this solution to the one obtained with the implementation of ridge regression from ``scikit-learn``:


```python
from sklearn.linear_model import Ridge
clf = Ridge(alpha=x.beta, fit_intercept=False)
y = np.ones(m)
clf.fit(X=A, y=y)
x_ridge = clf.coef_

difference_to_ridge = np.linalg.norm(x_min - x_ridge)
print(f"Difference to ridge: {difference_to_ridge}")
```

    Difference to ridge: 3.467185801042848e-14


The two solutions agree up to a precision of $10^{-13}$!

Note that the ``solution`` object also provides us access to the minimum of the cost function. Here, it is important to keep in mind that the cost function always has a factor $\frac{1}{2}$ in front.


```python
cost_at_minimum = solution.cost
print(f"Cost at minimum: {cost_at_minimum}")
```

    Cost at minimum: 100.0999000999005


Let us verify this by manually recomputing the cost at ``x_min``:


```python
cost_recomputed = 0.5 * np.sum(np.square(F(x_min))) + 0.5 * x.beta * np.sum(np.square(x_min))
print(f"Cost at minimum, recomputed: {cost_at_minimum}")
```

    Cost at minimum, recomputed: 100.0999000999005


And indeed, both numbers agree.

Finally, the ``solution`` also provides access to the [precision matrix](https://en.wikipedia.org/wiki/Precision_(statistics)) via ``solution.precision``.



```python
precision = solution.precision
precision.shape
```




    (200, 200)



This is relevant in the case where our optimization problem comes from [maximum likelihood estimation](https://en.wikipedia.org/wiki/Maximum_likelihood_estimation) or
[Bayesian maximum-a-posteriori estimation](https://en.wikipedia.org/wiki/Maximum_a_posteriori_estimation), and the cost function actually corresponds to a negative log-likelihood or a negative log-posterior density, respectively.

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