

<center>
    <h1>ILI285/INF285 Computación Científica </h1>
    <h1>Pauta Pregunta de "Newton + Gradiente Conjugado" - COP3</h1>
</center>

Necesitamos ajustar un conjunto de datos $D=\{(x_1,y_1),(x_2,y_2), \dots, (x_n, y_n)\}$, con una Spline cúbica de un intervalo, es decir, $S(x)=a+bx+cx^2+dx^3$. Para esto requerimos minimizar la función $F(a,b,c,d)=\displaystyle \sum_{i=1}^n (y_i- S(x_i))^4=\sum_{i=1}^n (y_i-a-bx_i-cx_i^2-dx_i^3)^4$. Esto significa que debemos obtener:
\begin{equation}
    \begin{split}
        \frac{\partial F}{\partial a} &= \sum_{i=1}^n -4(y_i-a-bx_i-cx_i^2-dx_i^3)^3 = 0 \\
        \frac{\partial F}{\partial b} &= \sum_{i=1}^n -4x_i(y_i-a-bx_i-cx_i^2-dx_i^3)^3 = 0 \\
        \frac{\partial F}{\partial c} &= \sum_{i=1}^n -4x_i^2(y_i-a-bx_i-cx_i^2-dx_i^3)^3 = 0 \\
        \frac{\partial F}{\partial d} &= \sum_{i=1}^n -4x_i^3(y_i-a-bx_i-cx_i^2-dx_i^3)^3 = 0 \\
    \end{split}
\end{equation}

Notar que simplificando la expresión anterior y definiendo $\mathbf{z}=(a,b,c,d)$ podemos llevar el resultado anterior al problema $\mathbf{G}(\mathbf{z})=\mathbf{0}$, donde $\mathbf{G}(\mathbf{z})$ se define como:

\begin{equation}
    \mathbf{G}(\mathbf{z}) = 
    \begin{bmatrix}
        \displaystyle \sum_{i=1}^n(y_i-a-bx_i-cx_i^2-dx_i^3)^3 \\
        \displaystyle \sum_{i=1}^nx_i(y_i-a-bx_i-cx_i^2-dx_i^3)^3 \\
        \displaystyle \sum_{i=1}^nx_i^2(y_i-a-bx_i-cx_i^2-dx_i^3)^3 \\
        \displaystyle \sum_{i=1}^nx_i^3(y_i-a-bx_i-cx_i^2-dx_i^3)^3
    \end{bmatrix} =
    \begin{bmatrix}
        0 \\ 0 \\ 0 \\ 0 
    \end{bmatrix}.
\end{equation}

Para obtener el $\mathbf{z}=(a,b,c,d)$ que minimiza $F$, modificaremos el método de Newton utilizando *Gradiente Conjugado* y así resolver el sistema de ecuaciones lineales que aparece en cada iteración. Si bien no sabemos de antemano si $J(\mathbf{z})$ es simétrica y definida positiva, $J(\mathbf{z})^TJ(\mathbf{z})$ si lo es, por lo que el nuevo sistema que se resuelve dentro del método de Newton es:
\begin{equation}
    J(\mathbf{z}_i)^TJ(\mathbf{z}_i) \Delta \mathbf{z}_i= -J(\mathbf{z}_i)^T\mathbf{G}(\mathbf{z}_i).
\end{equation}

¿Cuál es el valor del parámetro $p\in\{a, b, c, d\}$ luego de $k\in\{2,3,4\}$ iteraciones del método de Newton? 

Considere como *initial guess* el vector nulo tanto para el método del *Gradiente Conjugado* como para el *Método de Newton*. Además utilice como $10^{-10}$ como tolerancia para el método del *Gradiente Conjugado*.

# Solución propuesta


```python
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact
import ipywidgets as widgets
```

## Implementación Gradiente Conjugado


```python
def conjugateGradient(A, b, x_0=None, n=None, tol=1e-10):
    if n == None:
        n = b.shape[-1]
    X = np.zeros((n + 1, b.shape[0]))
    R = np.zeros_like(X)
    D = np.zeros_like(X)
    if x_0 is not None:
        X[0] = x_0
    R[0] = b - np.dot(A, X[0])
    D[0] = R[0]
    for k in range(n):
        a_k = np.dot(D[k], R[k]) / np.dot(D[k], np.dot(A, D[k]))
        X[k+1] = X[k] + a_k * D[k]
        R[k+1] = b - np.dot(A, X[k+1])
        b_k = np.dot(D[k], np.dot(A, R[k+1])) / np.dot(D[k], np.dot(A, D[k]))
        D[k+1] = R[k+1] - b_k * D[k]
        if np.linalg.norm(R[k+1]) < tol:
            X = X[:k+2]
            R = R[:k+2]
            D = D[:k+2]
            break
    return X[-1]#, R, D
```

## Implementación Newton en $\mathbb{R}^n$


```python
def newtonRn(F, J, x_0, n, tol=1e-10):
    x = np.zeros((n + 1, x_0.shape[0]))
    x[0] = x_0
    for k in range(n):
        JT = J(x[k]).T
        JTJ = np.dot(JT, J(x[k]))
        JTF = np.dot(JT, F(x[k]))
        w = conjugateGradient(JTJ, -JTF)
        x[k+1] = x[k] + w
        if np.linalg.norm(F(x[k+1])) < tol:
            x = x[:k+2]
            break
    return x
```

## Spline cúbica


```python
S = lambda a, b, c, d, x: a + b * x + c * x ** 2 + d * x ** 3 
```

## Cálculo del Jacobiano de $\mathbf{G}(\mathbf{z})$

La matriz Jacobiana asociada a este problema se define como:

\begin{equation}
    \scriptsize
    J(\mathbf{z})=
    \begin{bmatrix}
        \displaystyle -3\sum_{i=1}^n(y_i-a-bx_i-cx_i^2-dx_i^3)^2  & 
            \displaystyle -3\sum_{i=1}^nx_i(y_i-a-bx_i-cx_i^2-dx_i^3)^2 & 
            \displaystyle -3\sum_{i=1}^nx_i^2(y_i-a-bx_i-cx_i^2-dx_i^3)^2 & 
            \displaystyle -3\sum_{i=1}^nx_i^3(y_i-a-bx_i-cx_i^2-dx_i^3)^2 \\
        \displaystyle -3\sum_{i=1}^nx_i(y_i-a-bx_i-cx_i^2-dx_i^3)^2 & 
            \displaystyle -3\sum_{i=1}^nx_i^2(y_i-a-bx_i-cx_i^2-dx_i^3)^2 & 
            \displaystyle -3\sum_{i=1}^nx_i^3(y_i-a-bx_i-cx_i^2-dx_i^3)^2 & 
            \displaystyle -3\sum_{i=1}^nx_i^4(y_i-a-bx_i-cx_i^2-dx_i^3)^2 \\
        \displaystyle -3\sum_{i=1}^nx_i^2(y_i-a-bx_i-cx_i^2-dx_i^3)^2 & 
            \displaystyle -3\sum_{i=1}^nx_i^3(y_i-a-bx_i-cx_i^2-dx_i^3)^2 & 
            \displaystyle -3\sum_{i=1}^nx_i^4(y_i-a-bx_i-cx_i^2-dx_i^3)^2 & 
            \displaystyle -3\sum_{i=1}^nx_i^5(y_i-a-bx_i-cx_i^2-dx_i^3)^2 \\
        \displaystyle -3\sum_{i=1}^nx_i^3(y_i-a-bx_i-cx_i^2-dx_i^3)^2 & 
            \displaystyle -3\sum_{i=1}^nx_i^4(y_i-a-bx_i-cx_i^2-dx_i^3)^2 & 
            \displaystyle -3\sum_{i=1}^nx_i^5(y_i-a-bx_i-cx_i^2-dx_i^3)^2 & 
            \displaystyle -3\sum_{i=1}^nx_i^6(y_i-a-bx_i-cx_i^2-dx_i^3)^2 \\
    \end{bmatrix}
\end{equation}

Función para construir $\mathbf{G}(\mathbf{z})$ y $J(\mathbf{z})$.


```python
def buildGJ(x_i, y_i):
    g1 = lambda z: np.sum((y_i - S(*z, x_i)) ** 3)
    g2 = lambda z: np.sum(x_i * (y_i - S(*z, x_i)) ** 3)
    g3 = lambda z: np.sum(x_i ** 2 * (y_i - S(*z, x_i)) ** 3)
    g4 = lambda z: np.sum(x_i ** 3 * (y_i - S(*z, x_i)) ** 3)
    G = lambda z: np.array([g1(z), g2(z), g3(z), g4(z)])
    J = lambda z: -3 * np.array([
        [np.sum((y_i - S(*z, x_i)) ** 2), 
             np.sum(x_i * (y_i - S(*z, x_i)) ** 2), 
             np.sum(x_i ** 2 * (y_i - S(*z, x_i)) ** 2), 
             np.sum(x_i ** 3 * (y_i - S(*z, x_i)) ** 2)], 
        [np.sum(x_i * (y_i - S(*z, x_i)) ** 2), 
             np.sum(x_i ** 2 * (y_i - S(*z, x_i)) ** 2), 
             np.sum(x_i ** 3 * (y_i - S(*z, x_i)) ** 2), 
             np.sum(x_i ** 4 * (y_i - S(*z, x_i)) ** 2)], 
        [np.sum(x_i ** 2 * (y_i - S(*z, x_i)) ** 2), 
             np.sum(x_i ** 3 * (y_i - S(*z, x_i)) ** 2), 
             np.sum(x_i ** 4 * (y_i - S(*z, x_i)) ** 2), 
             np.sum(x_i ** 5 * (y_i - S(*z, x_i)) ** 2)], 
        [np.sum(x_i ** 3 * (y_i - S(*z, x_i)) ** 2), 
             np.sum(x_i ** 4 * (y_i - S(*z, x_i)) ** 2), 
             np.sum(x_i ** 5 * (y_i - S(*z, x_i)) ** 2), 
             np.sum(x_i ** 6 * (y_i - S(*z, x_i)) ** 2)]
    ])
    return G, J
```

# Respuestas

Función para obtener los parámetros de $S$ dado un dataset.


```python
def solution(x_i, y_i, k, GJ, n_par):
    x_0 = np.zeros(n_par)
    G, J = GJ(x_i, y_i)
    p = newtonRn(G, J, x_0, k, tol=1e-10)
    return p[-1]
```

### Lectura de datos


```python
DIR = 'data/'

# NPY
dataset_npy_1 = np.load(DIR + '1.npy')
dataset_npy_2 = np.load(DIR + '2.npy')
dataset_npy_3 = np.load(DIR + '3.npy')
y_npy_1 = dataset_npy_1[:, 1]
y_npy_2 = dataset_npy_2[:, 1]
y_npy_3 = dataset_npy_3[:, 1]

# CSV
dataset_csv_1 = np.loadtxt(DIR + '1.csv', delimiter=',')
dataset_csv_2 = np.loadtxt(DIR + '2.csv', delimiter=',')
dataset_csv_3 = np.loadtxt(DIR + '3.csv', delimiter=',')
y_csv_1 = dataset_csv_1[:, 1]
y_csv_2 = dataset_csv_2[:, 1]
y_csv_3 = dataset_csv_3[:, 1]

# TXT
dataset_txt_1 = np.loadtxt(DIR + '1.txt', delimiter=',')
dataset_txt_2 = np.loadtxt(DIR + '2.txt', delimiter=',')
dataset_txt_3 = np.loadtxt(DIR + '3.txt', delimiter=',')
y_txt_1 = dataset_txt_1[:, 1]
y_txt_2 = dataset_txt_2[:, 1]
y_txt_3 = dataset_txt_3[:, 1]
```

Selección de datos. Salvo el archivo ```2.npy```, todos los otros dataset son iguales solo cambia el formato del archivo.


```python
n = 100
x_a, x_b = -1, 1
x_i = np.linspace(x_a, x_b, n)
y_i1 = dataset_npy_1[:, 1] # 1.{npy, csv, txt}
y_i2 = dataset_npy_2[:, 1] # 2.npy
y_i22 = dataset_csv_2[:, 1] # 2.{csv, txt}
y_i3 = dataset_npy_3[:, 1] # 3.{npy, csv, txt}
```

Combinaciones de parámetros.


```python
K = [2, 3, 4] # Number of iterations
IP = [0, 1, 2, 3] # Parameter index
PL = ["a", "b", "c", "d"] # Parameter name
D = [y_i1, y_i2, y_i22, y_i3] # Dataset
Dn = ['1.{npy, csv, txt}', '2.npy', '2.{csv, txt}', '3.{npy, csv, txt}'] # Dataset name
```


```python
def experiment(i, p, k):
    par = solution(x_i, D[i], k, buildGJ, 4)
    print("k: %d, parametro: %s, valor: %f" % (k, PL[p], par[p]))
    
interact(experiment,
        i=widgets.Dropdown(
            options=[('1.{npy, csv, txt}', 0), ('2.npy', 1), ('2.{csv, txt}', 2), ('3.{npy, csv, txt}', 3)],
            value=0,
            description='Dataset:'
        ),
        p=widgets.Dropdown(
            options=[("a", 0), ("b", 1), ("c", 2), ("d", 3)],
            value=0,
            description='Parámetro:'
        ),
        k=K
)
```


    interactive(children=(Dropdown(description='Dataset:', options=(('1.{npy, csv, txt}', 0), ('2.npy', 1), ('2.{c…





    <function __main__.experiment(i, p, k)>



5ed21573156c96503617e27c612426c67c638652
ipynb
Jupyter Notebook
cop-3/newton_gc.ipynb
fitoxgs2/INF-285
f7d547c2480aea7eb1242b873ab1135a0873834f
BSD-3-Clause
cop-3/newton_gc.ipynb
fitoxgs2/INF-285
f7d547c2480aea7eb1242b873ab1135a0873834f
BSD-3-Clause
cop-3/newton_gc.ipynb
fitoxgs2/INF-285
f7d547c2480aea7eb1242b873ab1135a0873834f
BSD-3-Clause
Qwen/Qwen-72B
1. YES
2. YES