```python
import numpy as np
import numpy.linalg as la
import sympy as sp
```


```python
def gradient(formula, symbols, values=None):
    '''
    Given a SymPy formula and variables
    Find its analytic gradient without substituion
    as a list of SymPy formulae or numerical gradient
    if values specified
    '''
    size = len(symbols)
    gradient = []
    for i in range(size):
        gradient.append(formula.diff(symbols[i]))
        
    if values == None:
        return gradient

    # Make sure you don't mess up the analytical gradient
    g_copy = gradient.copy()
    gradient_at = []
    for i in range(len(g_copy)):
        for j in range(len(symbols)):
            g_copy[i] = g_copy[i].subs(symbols[j], values[j])
        gradient_at.append(float(g_copy[i].evalf()))
    return gradient_at
```


```python
def subs_all(formula, variables, values):
    '''
    You know what, it's getting to the point where this function
    is necessary
    '''
    result = formula
    for i in range(len(values)):
        result = result.subs(variables[i], values[i])
    return float(result.evalf())
def hessian(gradient, symbols, values=None):
    '''
    Given an analytic gradient and variables
    Calculate its analytic Hessian
    or numerical Hessian if values are specified
    '''
    size = len(symbols)
    hessian = []
    for i in range(size):
        hessian.append([0]*size)
    for i in range(size):
        for j in range(size):
            hessian[i][j] = gradient[i].diff(symbols[j])
            
    if values == None:
        return hessian
    
    hessian_at = hessian.copy()
    size = len(hessian)
    for i in range(size):
        for j in range(size):
            hessian_at[i][j] = subs_all(hessian_at[i][j], symbols, [float(v) for v in values])
#             for k in range(len(symbols)):
#                 hessian_at[i][j] = hessian_at[i][j].subs(symbols[k], values[k])
#             hessian_at[i][j] = float(hessian_at[i][j].evalf())
    return hessian_at
```


```python
def newton_nd_optimization_crude_step(formula, symbols, x_prev):
    x_prev = list(x_prev)
    grad = gradient(formula, symbols)
    neg_gradient_at = -1 * np.array(gradient(formula, symbols, values=x_prev))
    hes_at = np.array(hessian(grad, symbols, values=x_prev))
    return np.array(x_prev) + la.solve(hes_at, neg_gradient_at)
```


```python
def newton_nd_optimization_crude(f_str, s_str, start, tolerance, actual_solution):
    '''
    A crude version of Newton ND Optimization algorithm
    Maybe you can help out implementing an analytic solution
    Returning the numerical solution as well as the number of
    iterations it took
    '''
    formula = sp.sympify(f_str)
    symbols = sp.symbols(s_str)
    curr = np.copy(start)
    iterations = 0
    while (la.norm(curr - actual_solution, 2) > tolerance):
        curr = newton_nd_optimization_crude_step(formula, symbols, curr)
        iterations += 1
    return curr, iterations
```


```python
# Test case set-up
form = sp.sympify('12*x**2+10*x*y+12*y**2+8*E**(9*x*y)+8*(sin(y)**2)+9*cos(x*y)') # 记住，大E才是natural number
x, y = sp.symbols('x y')
tolerance = 10**-7
# Test Case: Gradient
grad_at = gradient(form, (x, y), values=[0, 1])
expected = np.array([
    float(form.diff(x).subs(x, 0).subs(y, 1).evalf()),
    float(form.diff(y).subs(x, 0).subs(y, 1).evalf())
])
actual = np.array(grad_at)
assert la.norm(expected - actual, 2) < tolerance
```


```python
# Test Case: Hessian
hes = hessian(gradient(form, (x, y)), (x, y), values=[0, 1])
expected = np.array([
    [
        float(form.diff(x).diff(x).subs(x, 0).subs(y, 1).evalf()),
        float(form.diff(x).diff(y).subs(x, 0).subs(y, 1).evalf())
    ], [
        float(form.diff(y).diff(x).subs(x, 0).subs(y, 1).evalf()),
        float(form.diff(y).diff(y).subs(x, 0).subs(y, 1).evalf())
    ]
])
actual = np.array(hes)
assert la.norm(expected - actual, 2) < tolerance
```


```python
# Test Case: Newton ND step
H = np.array([
    [663, 82],
    [82, -16*np.sin(1)**2 + 16*np.cos(1)**2 + 24]
])
v = np.array([82, 16*np.sin(1)*np.cos(1) + 24])
expected = la.solve(H, v) * -1 + np.array([0, 1])
actual = newton_nd_optimization_crude_step(form, (x, y), np.array([0, 1]))
assert la.norm(expected - actual, 2) < tolerance
```


```python
form = sp.sympify('11*x**2+13*x*y+11*y**2+6*E**(2*x*y)+5*(sin(y)**2)+7*cos(x*y)')
symbols = sp.symbols('x y')
print(newton_nd_optimization_crude_step(form, symbols, [0, 1]))
print(gradient(form, symbols, values=[0, 1]))
print(hessian(gradient(form, symbols), symbols, values=[0, 1]))
# Works, great!
```

    [ 3.07907313 -4.80335408]
    [25.0, 26.54648713412841]
    [[39.0, 25.0], [25.0, 17.83853163452858]]



```python
# Workspace
form = sp.sympify('4*x**2+6*y**4')
symbols = sp.symbols('x y')
values = [2, 1]
hes = hessian(gradient(form, symbols), symbols, values)
x_1 = newton_nd_optimization_crude_step(form, symbols, values)
step = np.array(x_1) - np.array(values)
print(hes)
print(x_1)
print(step)
```

    [[8.0, 0.0], [0.0, 72.0]]
    [0.         0.66666667]
    [-2.         -0.33333333]



```python

```

b80fc7e2eb5099b521be47ad0edb3509b97baa4e
ipynb
Jupyter Notebook
newton_nd_optimization_crude.ipynb
Racso-3141/uiuc-cs357-fa21-scripts
e44f0a1ea4eb657cb77253f1db464d52961bbe5e
MIT
2021-11-02T05:56:10.000Z
2022-03-03T19:25:19.000Z
newton_nd_optimization_crude.ipynb
Racso-3141/uiuc-cs357-fa21-scripts
e44f0a1ea4eb657cb77253f1db464d52961bbe5e
MIT
newton_nd_optimization_crude.ipynb
Racso-3141/uiuc-cs357-fa21-scripts
e44f0a1ea4eb657cb77253f1db464d52961bbe5e
MIT
2021-10-30T15:18:01.000Z
2021-12-10T11:26:43.000Z
Qwen/Qwen-72B
1. YES
2. YES