

<center>
    <h1>ILI285/INF285 Computación Científica </h1>
    <h1> COP-4</h1>
    <h1> LM: Levenberg-Marquardt</h1>
</center>

# Pregunta

Considere el siguiente conjunto de datos $D=\{(x_1,y_1),(x_2,y_2), \dots, (x_n, y_n)\}$ que se quieren ajustar con la siguiente función $f(x, \boldsymbol{\beta})=\beta_1+\sin(\beta_2\,x)+\cos(\beta_3\,x)$, donde $\boldsymbol{\beta}=(\beta_1, \beta_2, \beta_3 )$.  Para esto requerimos minimizar la función $S(\boldsymbol{\beta})=\displaystyle \sum_{i=1}^n \left(y_i- f(x_i, \boldsymbol{\beta})\right)^{\{2, 4\}}$. De similar forma a lo realizado en el certamen anterior, debemos encontrar un mínimo, para lo cual debemos obtener las derivadas parciales con respecto a los parámetros, esto es:

\begin{equation}
    \begin{split}
        \frac{\partial S}{\partial \beta_1} &= 
        \sum_{i=1}^n -\{2, 4\}(y_i-\beta_1-\sin(\beta_2\,x_i)-\cos(\beta_3\,x_i))^{\{1, 3\}} = 0 \\
        \frac{\partial S}{\partial \beta_2} &= 
        \sum_{i=1}^n -\{2, 4\}x_i\cos(\beta_2\,x_i)(y_i-\beta_1-\sin(\beta_2\,x_i)-\cos(\beta_3\,x_i))^{\{1, 3\}} = 0 \\
        \frac{\partial S}{\partial \beta_3} &= 
        \sum_{i=1}^n \{2, 4\}x_i\sin(\beta_3\,x_i)(y_i-\beta_1-\sin(\beta_2\,x_i)-\cos(\beta_3\,x_i))^{\{1, 3\}} = 0 
    \end{split}
\end{equation}


Podemos llevar la expresión anterior al problema $\mathbf{F}(\boldsymbol{\beta})=\mathbf{0}$, donde $\mathbf{F}(\boldsymbol{\beta})$ se define como:

\begin{equation}
    \mathbf{F}(\boldsymbol{\beta}) = 
    \begin{bmatrix}
        \displaystyle \sum_{i=1}^n(y_i-\beta_1-\sin(\beta_2\,x_i)-\cos(\beta_3\,x_i))^{\{1, 3\}} \\
        \displaystyle \sum_{i=1}^nx_i\cos(\beta_2\,x_i)\,(y_i-\beta_1-\sin(\beta_2\,x_i)-\cos(\beta_3\,x_i))^{\{1, 3\}} \\
        \displaystyle \sum_{i=1}^nx_i\sin(\beta_3\,x_i)\,(y_i-\beta_1-\sin(\beta_2\,x_i)-\cos(\beta_3\,x_i))^{\{1, 3\}}
    \end{bmatrix}. \tag{1}
\end{equation}


Para obtener el $\boldsymbol{\beta}$ que minimiza $S$, utilizaremos una variante del método de Newton en alta dimensión, el cual es conocido como el método de Levenberg-Marquardt. En resumen, este método es muy similar al método de Newton en $\mathbb{R}^n$ pero el sistema de ecuaciones que se resuelve en cada iteración es el siguiente:

\begin{equation}
    \begin{split}
        \left(J(\boldsymbol{\beta}_i)^TJ(\boldsymbol{\beta}_i) + 
        \lambda\,I\right) \Delta \boldsymbol{\beta}_i
        & = -J(\boldsymbol{\beta}_i)^T\mathbf{F}(\boldsymbol{\beta}_i),\\
        \boldsymbol{\beta}_{i+1}&=\boldsymbol{\beta}_i+\Delta \boldsymbol{\beta}_i
    \end{split}\tag{2}
\end{equation}


donde $J(\boldsymbol{\beta}_i)$ es la matriz Jacobiana de $\mathbf{F}$ evaluada en $\boldsymbol{\beta}_i$. El objetivo del parámetro $\lambda$ es reducir el número de condición asociado a la matriz $J(\boldsymbol{\beta}_i)^T\,J(\boldsymbol{\beta}_i)$, además de mejorar la convergencia del método. Sin embargo, al utilizar un valor de $\lambda$ muy grande, la solución obtenida es una perturbación de la solución límite, es decir, cuando $\lambda$ tiende a $0$.  Para evitar calcular la matriz Jacobiana $J(\boldsymbol{\beta})$ utilizaremos una aproximación. Considere la siguiente linealización de la función $\mathbf{F}$:

\begin{equation}
\mathbf{F}(\boldsymbol{\beta} + \varepsilon\,\mathbf{v}) = \mathbf{F}(\boldsymbol{\beta}) + 
J(\boldsymbol{\beta})(\varepsilon\,\mathbf{v}) + \text{Términos de orden } \varepsilon^2 + \dots.
\end{equation}


donde $\mathbf{v}$ es un vector que será definido más adelante y $0 <\varepsilon \ll 1$ un escalar.  Basado en la expansión anterior, podemos despejar el producto de la matriz Jacobiana por el vector $\mathbf{v}$ de la siguiente forma:

\begin{equation}
    J(\boldsymbol{\beta})\,\mathbf{v} \approx
    \frac{\mathbf{F}(\boldsymbol{\beta} + \varepsilon\,\mathbf{v}) 
    - \mathbf{F}(\boldsymbol{\beta})}{\varepsilon}\tag{3}
\end{equation}


Con esta aproximación, podemos notar que al definir $\mathbf{v}=\mathbf{e}_1=[1,0,0]^T$, se obtiene una aproximación de la primera columna de $J(\boldsymbol{\beta})$, con $\mathbf{v}=\mathbf{e}_2=[0,1,0]^T$ la segunda columna y $\mathbf{v}=\mathbf{e}_3=[0,0,1]^T$ la tercera columna. Dicho de otra manera, tenemos:

\begin{equation}
    J(\boldsymbol{\beta}) \approx
    \left[
        \begin{array}{c|c|c}
            & & \\
            \frac{\mathbf{F}(\boldsymbol{\beta} + \varepsilon\,\mathbf{e}_1) 
            - \mathbf{F}(\boldsymbol{\beta})}{\varepsilon} & 
            \frac{\mathbf{F}(\boldsymbol{\beta} + \varepsilon\,\mathbf{e}_2) 
            - \mathbf{F}(\boldsymbol{\beta})}{\varepsilon} & 
            \frac{\mathbf{F}(\boldsymbol{\beta} + \varepsilon\,\mathbf{e}_3) 
            - \mathbf{F}(\boldsymbol{\beta})}{\varepsilon} \\
            & &
        \end{array}
    \right]
\end{equation}

 

La gran ventaja de la aproximación anterior es que uno puede obtener numéricamente una estimación razonable de la matriz Jacobiana omitiendo el manejo algebraico de obtener el gradiente de cada termino en la ecuación $(1)$.

Considerando el cálculo de la matriz Jacobiana aproximada por la ecuación $(3)$ con $\varepsilon=10^{-10}$ y la iteración de Levenberg-Marquardt definida en $(2)$, determine el valor del parámetro $\beta_1$ luego de $\{3,4,5\}$ iteraciones del método de Levenberg-Marquardt utilizando $\lambda = 1$ y para resolver el sistema de ecuaciones en $(2)$, utilice al algoritmo de {Jacobi, Gauss-Seidel, SOR($1.1$)} con $3$ iteraciones, *initial guess* $\boldsymbol{\beta}_0=[1,1,1]^T$ y $10^{-10}$ como tolerancia en {Jacobi, Gauss-Seidel, SOR($1.1$)}.

Los datos para ajustar el modelo se encuentran en https://github.com/sct-utfsm/INF-285/tree/master/cop/4/data/lm/1.npy. En la primera columna encontrará los valores de $x_i$ y la segunda columna los valores de $y_i$.

**Debe entregar sus resultados con 5 decimales sin redondear. Por ejemplo si su resultado es 2365.1345871237 debe completar con 2365.13458.**

**Preguntas:**

1. ¿Cuál es el valor de la segunda componente de $\mathbf{F}(\boldsymbol{\beta}_0)$? 
2. ¿Cuál es la norma de Frobenius de la aproximación del Jacobiano en $\boldsymbol{\beta}_0$, es decir,  $\|J(\boldsymbol{\beta}_0)\|_F$? 
3. ¿Cuál es el valor del parámetro $\beta_1$ luego de ajustar el modelo? 

## Bibliotecas


```python
import numpy as np
import scipy.linalg as spla
import matplotlib.pyplot as plt
```

## *Solvers* de sistemas de ecuaciones lineales


```python
def jacobi(A, b, x_0=None, max_iter=50, tol=1e-10):
    n = A.shape[0] # Matrix size
    X = np.zeros((max_iter + 1, n)) # Matrix with solution at each iteration
    # Initial guess
    if x_0 is not None:
        X[0] = x_0
    D = np.diag(A) # Diagonal of A (only keep a vector with diagonal)
    D_inv = np.diag(1 / D) # Inverse of D
    r = b - np.dot(A, X[0]) # Residual vector
    # Jacobi iteration
    for k in range(max_iter):
        X[k+1] = X[k] + np.dot(D_inv, r)
        r = b - np.dot(A, X[k+1]) # Update residual
        if np.linalg.norm(r) < tol: # Stop criteria
            X = X[:k+2]
            break
    return X[-1] 
```


```python
def gaussSeidel(A, b, x_0=None, max_iter=50, tol=1e-10):
    n = A.shape[0] # Matrix size
    X = np.zeros((max_iter + 1, n)) # Matrix with solution at each iteration
    # Initial guess
    if x_0 is not None:
        X[0] = x_0
    LD = np.tril(A) # Get lower triangle (L + D)
    # Get inverse in O(n^2) instead of np.linalg.inv(LD) O(n^3)
    LD_inv = spla.solve_triangular(LD, np.eye(n), lower=True) 
    r = b - np.dot(A, X[0]) # Residual
    # Gauss-Seidel iteration
    for k in range(max_iter):
        X[k+1] = X[k] + np.dot(LD_inv, r)
        r = b - np.dot(A, X[k+1]) # Residual update
        if np.linalg.norm(r) < tol: # Stop criteria
            X = X[:k+2]
            break
    return X[-1]
```


```python
def SOR(A, b, x_0=None, w=1.05, max_iter=50, tol=1e-10):
    n = A.shape[0] # Matrix size
    X = np.zeros((max_iter + 1, n)) # Matrix with solution at each iteration
    # Initial guess
    if x_0 is not None:
        X[0] = x_0
    L = np.tril(A, k=-1) # Get lower triangle 
    Dw = np.diag(np.diag(A) / w)
    # Get inverse in O(n^2) instead of np.linalg.inv(L+Dw) O(n^3)
    LDw_inv = spla.solve_triangular(L + Dw, np.eye(n), lower=True) 
    r = b - np.dot(A, X[0]) # Residual
    # SOR iteration
    for k in range(max_iter):
        X[k+1] = X[k] + np.dot(LDw_inv, r)
        r = b - np.dot(A, X[k+1]) # Residual update
        if np.linalg.norm(r) < tol: # Stop criteria
            X = X[:k+2]
            break
    return X[-1]
```

## Levenberg-Marquardt


```python
def levenbergMarquardt(F, J, x_0, l, solver, max_iter=100, tol=1e-10):
    x = np.zeros((max_iter + 1, x_0.shape[0]))
    x[0] = x_0
    for k in range(max_iter):
        JT = J(x[k]).T # Compute transpose of J
        JTJ = np.dot(JT, J(x[k]))  # Compute J.T J
        JTJ_l = JTJ + l * np.eye(x_0.shape[0])  # J.T J + \lambda I
        JTF = np.dot(JT, F(x[k])) 
        w = solver(JTJ_l, -JTF)
        x[k+1] = x[k] + w
        if np.linalg.norm(F(x[k+1])) < tol:
            x = x[:k+2]
            break
    return x
```


```python
f = lambda x, beta: beta[0] + np.sin(beta[1] * x) + np.cos(beta[2] * x) 
```


```python
def FD(F, x, v, eps=1e-10):
    return (F(x + eps * v) - F(x)) / eps
```


```python
def buildFJ(x_i, y_i, f_opt):
    f1 = f_opt['f1']
    f2 = f_opt['f2']
    f3 = f_opt['f3']
    F = lambda beta: np.array([f1(beta), f2(beta), f3(beta)])
    I = np.eye(3)
    J = lambda beta: np.array([
        FD(F, beta, I[0]), FD(F, beta, I[1]), FD(F,beta, I[2])
    ]).T
    return F, J
```


```python
data = np.load('data/lm/1.npy')
x_i, y_i = data[:, 0], data[:, 1]
```


```python
def experiment(Sopt, sol, lam, solv_iter, lm_iter, lm_b0, sol_x0, plo, verb=0):
    if lm_b0 == sol_x0 == 0:
        print("Opcion no valida")
        return None
    
    b0 = np.ones(3) * lm_b0
    x0 = np.ones(3) * sol_x0
    if Sopt == 0:
        f_S = {
            'f1': lambda beta: np.sum(y_i - f(x_i, beta)),
            'f2': lambda beta: np.sum(x_i * np.cos(beta[1] * x_i) * (y_i - f(x_i, beta))),
            'f3': lambda beta: np.sum(x_i * np.sin(beta[2] * x_i) * (y_i - f(x_i, beta)))
        }
    elif Sopt == 1:
        f_S = {
            'f1': lambda beta: np.sum((y_i - f(x_i, beta)) ** 3),
            'f2': lambda beta: np.sum(x_i * np.cos(beta[1] * x_i) * (y_i - f(x_i, beta)) ** 3),
            'f3': lambda beta: np.sum(x_i * np.sin(beta[2] * x_i) * (y_i - f(x_i, beta)) ** 3)
        }
        
    F, J = buildFJ(x_i, y_i, f_S)
    
    F_b0 = F(b0)
    J_b0 = J(b0)
    
    if verb:
        print("Respuestas: ")
        print("--------------------")
        print("1. Segunda componente de F(\\vec{beta}_0): %.6f" %  F_b0[1])
        print("2. ||J(\\vec{beta}_0)||_F: %.6f" % np.linalg.norm(J_b0, ord="fro"))

    ss = [jacobi, gaussSeidel, SOR]
    
    if sol < 2:
        solver = lambda A, b: ss[sol](A, b, max_iter=solv_iter, x_0=x0) 
    elif sol == 2:
        solver = lambda A, b: ss[sol](A, b, max_iter=solv_iter, w=1.1, x_0=x0) 
        
    B = levenbergMarquardt(F, J, b0, lam, solver, lm_iter, tol=1e-10)
    
    if verb:
        print("3. beta_1 luego de %d iteraciones: %.6f" %(lm_iter, B[-1, 0]))
    if plo:
        plt.plot(x_i, y_i, 'r.')
        plt.plot(x_i, f(x_i, B[-1]), '--', label='Curve')
        plt.grid(True)
        plt.xlabel(r'$x$')
        plt.ylabel(r'$y$')
        plt.legend()
        plt.show()
```

# Respuestas

Se considera el uso de:
* $\boldsymbol{\beta}_0=[1,1,1]$ como *initial guess* tanto para Levenberg-Marquard como para el *Solver* iterativo.
* $\boldsymbol{\beta}_0=[1,1,1]$ como *initial guess* de Levenberg-Marquard y $[0,0,0]$ como *initial guess*  para el *Solver* iterativo.
* $\boldsymbol{\beta}_0=[0,0,0]$ como *initial guess* de Levenberg-Marquard y $[1,1,1]$ como *initial guess*  para el *Solver* iterativo.


```python
from ipywidgets import *
```


```python
# widgets
sol = Dropdown(options=[('Jacobi', 0), ('Gauss-Seidel', 1), ('SOR(1.1)', 2)], value=0, description="Solver")
Sopt = Dropdown(options=[('S1', 0), ('S2', 1)], value=0, description="Función de Error")
lit = Dropdown(options=[3, 4, 5], description="Levenberg-Marquardt solver")
lm_b0_ = Dropdown(options=[('[0, 0, 0]', 0), ('[1, 1, 1]', 1)], 
                  value=1, description=r"L-M: $\boldsymbol{\beta}_0$")
sol_x0_ = Dropdown(options=[('[0, 0, 0]', 0), ('[1, 1, 1]', 1)], 
                   value=1, description=r"Solver: $\boldsymbol{\beta}_0$")        
plo = Checkbox(
    value=False,
    description='Plot',
    disabled=False,
    indent=False
)
```


```python
def onChange1(args):
    if args['new'] == 0:
        sol_x0_.value = 1
        
def onChange2(args):
    if args['new'] == 0:
        lm_b0_.value = 1
        
lm_b0_.observe(onChange1, 'value')
sol_x0_.observe(onChange2, 'value')

interact(experiment, Sopt=Sopt, sol=sol, lam=fixed(1), solv_iter=fixed(3), 
         lm_iter=lit, lm_b0=lm_b0_, sol_x0=sol_x0_, plo=plo, verb=fixed(1))
```


    interactive(children=(Dropdown(description='Función de Error', options=(('S1', 0), ('S2', 1)), value=0), Dropd…





    <function __main__.experiment(Sopt, sol, lam, solv_iter, lm_iter, lm_b0, sol_x0, plo, verb=0)>



036527d7891bb13608558dc080c8ae9e35f2c0c6
ipynb
Jupyter Notebook
cop/4/lm.ipynb
DoVirgo/INF-285
a2d994e713c89052820919aeb5ae4d39546599c0
BSD-3-Clause
2020-04-24T01:25:24.000Z
2021-05-10T01:08:37.000Z
cop/4/lm.ipynb
DoVirgo/INF-285
a2d994e713c89052820919aeb5ae4d39546599c0
BSD-3-Clause
2020-04-22T00:57:29.000Z
2020-04-23T23:59:28.000Z
cop/4/lm.ipynb
DoVirgo/INF-285
a2d994e713c89052820919aeb5ae4d39546599c0
BSD-3-Clause
2020-04-20T01:15:09.000Z
2020-07-27T22:53:59.000Z
Qwen/Qwen-72B
1. YES
2. YES