# Multivariate function and its derivative


```python
import torch
from torch.autograd import Variable
from torch.autograd import Function
import torch.nn.functional as F
import numpy as np

x = Variable(torch.tensor([3.0,-1.0,0.0,1.0]),requires_grad=True)
f = (x[0]+10.0*x[1])**2 + 5.0*(x[2]-x[3])**2 + (x[1]+2.0*x[2])**4 + 10.0*(x[0]-x[3])**4

print("Vector variable x:",x.data)

print("Function f at x:",f.data)

# compute gradient of f at x
g = torch.autograd.grad(f,x)

print("Gradient of f at x:",g[0].data)
```

    Vector variable x: tensor([ 3., -1.,  0.,  1.])
    Function f at x: tensor(215.)
    Gradient of f at x: tensor([ 306., -144.,  -18., -310.])


### Let's verify with math formula.

$f(x_1,x_2,x_3,x_4) = (x_1+10x_2)^2 + 5(x_3-x_4)^2 + (x_2+2x_3)^4 + 10(x_1-x_4)^4$

Let's evaluate these partial derivatives at

$x_1 = 3, x_2 = -1, x_3 = 0, x_4 = 1$

$\frac{\partial f}{\partial x_1} = 2(x_1+10x_2)+40(x_1-x_4)^3 = 2(3-10)+40(3-1)^3 = -14+320 = 306$

$\frac{\partial f}{\partial x_2} = 20(x_1+10x_2)+4(x_2+2x_3)^3 = 20(3-10)+4(-1)^3 = -140-4 = -144$

$\frac{\partial f}{\partial x_3} = 10(x_3-x_4)+8(x_2+2x_3)^3 = 10(-1)+8(-1)^3 = -10-8 = -18$

$\frac{\partial f}{\partial x_4} = -10(x_3-x_4)-40(x_1-x_4)^3 = -10(-1)-40(3-2)^3 = 10-320 = -310$

So, we have

$\nabla f(3,-1,0,1) = \left[ 306,-144,-18,-310 \right]$


# Numerical Derivative
### For any function we can also compute its numerical derivative, which is an approximation. The formula for computing numerical derivatives often uses Taylor series approximation from Calculus:
$\frac{\partial f}{\partial x_1} = \frac{f(x_1+\delta x_1,x_2,x_3,...)-f(x_1-\delta x_1,x_2,x_3,...)}{2\delta x_1},$

$\frac{\partial f}{\partial x_2} = \frac{f(x_1,x_2+\delta x_2,x_3,...)-f(x_1,x_2-\delta x_2,x_3,...)}{2\delta x_2},$

$...,$

### and so on for each variable $x_i$. As a result, comoutation of numerical derivative is extremely inefficient for function with a large number of variables. But it is often used as a check to see if your autograd code is working correctly, or not.



```python
x = Variable(torch.tensor([3.0,-1.0,0.0,1.0]),requires_grad=True)

def compute_f(x):
  f = (x[0]+10.0*x[1])**2 + 5.0*(x[2]-x[3])**2 + (x[1]+2.0*x[2])**4 + 10.0*(x[0]-x[3])**4
  return f

f = compute_f(x)

print("Vector variable x:",x.data)

print("Function f at x:",f.data)

# compute gradient of f at x
g = torch.autograd.grad(f,x)

print("Exact gradient of f at x:",g[0].data)

num_g = torch.zeros(4)
h=1e-3
eye = torch.eye(4) # 4-by-4 identity matrix
for i in range(4):
  num_g[i] = compute_f(x+h*eye[i,:]) - compute_f(x-h*eye[i,:])

num_g = num_g/(2.0*h)

print("Numerical gradient of f at x:",num_g.data)
```

    Vector variable x: tensor([ 3., -1.,  0.,  1.])
    Function f at x: tensor(215.)
    Exact gradient of f at x: tensor([ 306., -144.,  -18., -310.])
    Numerical gradient of f at x: tensor([ 305.9769, -143.9972,  -18.0054, -309.9976])


# PyTorch autograd has no problem working with branches!


```python
def compute_func(x):
  if x[0]<=0 and x[1]>=2:
    return x[0]**2 + x[1]**2Positive wh
  else:
    return -2.0*x[0]

x = Variable(torch.tensor([-1.0,3.0]),requires_grad=True)

f = compute_func(x)
g = torch.autograd.grad(f,x)
print(g)
```

    (tensor([-2.,  6.]),)


# Autograd can work with loops too!

### Here is an example. The function (due to Kepler https://en.wikipedia.org/wiki/Kepler%27s_equation#Inverse_Kepler_equation) below

$y = f(x) = x - 0.2sin(x)$

### has no closed form solution for its inverse $x = f^{-1}(y)$

### We can compute an approximate inverse $x = g(y)$ by an iterative method

### From calculus we know that $\frac{df^{-1}}{dy}= \frac{1}{\frac{df}{dx}}$

### So, derivative of $g$ computed by autograd can be verifed against calculus computed derivative: $\frac{dg}{dy}\approx \frac{1}{\frac{df}{dx}}$.




```python
def compute_Kepler(x):
    y = x - 0.9*torch.sin(x)
    return y

# Apply Newton's method to compute inverse
def compute_invKepler(y):
    x = y.detach().clone()
    for i in range(100):
        #x = x - (x-0.9*torch.sin(x)-y)/(1-0.9*torch.cos(y))
        x = y + 0.9*torch.sin(x) # fixed point iteration
    return x

x = Variable(torch.tensor([20.]),requires_grad=True)
y = compute_Kepler(x)
print("x:",x.item(),"y:",y.item())

x_ = compute_invKepler(y)

print("y:",y.item(),"x_:",x_.item())
print("So, inverse apprximation is fine.")

# compute derivative of x at y using autograd on inverse
dx_ = torch.autograd.grad(x_,y) # Note that this autograd has to take care of the for loop
print("dx_/dy by Autograd that needs to go through the for loop:",dx_[0])

# compute derivative of x at y
dy_dx = torch.autograd.grad(y,x) # this can also be computed using math formula: f'(x) = 1 - 0.5cos(x)
print("dx/dy by Calculus:",1./dy_dx[0])


```

    x: 20.0 y: 19.178348541259766
    y: 19.178348541259766 x_: 19.999998092651367
    So, inverse apprximation is fine.
    dx_/dy by Autograd that needs to go through the for loop: tensor([1.5805])
    dx/dy by Calculus: tensor([1.5805])


# Now that we know Calculs works(!), can we write our own derivative instead of using autograd?

## By the way, what might be an advantage of writing your own derivative here?


```python
# Inherit from Function
class My_invKepler(Function):

    # Note that both forward and backward are @staticmethods
    @staticmethod
    def forward(ctx, y):
        x = y
        for i in range(50):
            #x = x - (x-0.9*torch.sin(x)-y)/(1-0.9*torch.cos(y))
            x = y + 0.9*torch.sin(x)
        # Save for derivative computation
        ctx.save_for_backward(x,y)
        return x

    @staticmethod
    def backward(ctx, output_grad):
        # retrieve saved tensors and use them in derivative calculation
        x, y = ctx.saved_tensors

        # Return Jacobian-vector product (chain rule) - more to say about chain later!
        input_grad = (1.0/(1.0-0.9*torch.cos(x)))*output_grad
        #print("Output grad",output_grad)
        return input_grad
    
x__ = My_invKepler.apply(y)
print("y:",y.data,"x__:",x__.data)

dx__ = torch.autograd.grad(x__,y)
print(dx__)
```

    y: tensor([19.1783]) x__: tensor([20.0000])
    (tensor([1.5805]),)


# Gradient descent


```python
x = Variable(torch.tensor([3.0,-1.0,0.0,1.0]),requires_grad=True)
  
steplength = 1e-3 # for gradient descent
for i in range(1000):
  # function
  f = (x[0]+10.0*x[1])**2 + 5.0*(x[2]-x[3])**2 + (x[1]+2.0*x[2])**4 + 10.0*(x[0]-x[3])**4
  # compute grdaient
  g = torch.autograd.grad(f,x)
  # adjust variable
  x = x - steplength*g[0]
  if i%100==0:
    print("Current variable value:",x.detach().numpy(),"Current function value:", f.item())

```

    Current variable value: [ 2.694 -0.856  0.018  1.31 ] Current function value: 215.0
    Current variable value: [ 1.7774506  -0.19877124  0.5913045   1.3086869 ] Current function value: 4.060643672943115
    Current variable value: [ 1.4662572  -0.16261245  0.5297499   1.0466366 ] Current function value: 2.330897808074951
    Current variable value: [ 1.2376873  -0.13571279  0.47490963  0.85692835] Current function value: 1.4002106189727783
    Current variable value: [ 1.0653068  -0.11562601  0.4296491   0.7175411 ] Current function value: 0.8787093758583069
    Current variable value: [ 0.9327728  -0.10033523  0.39197212  0.6131387 ] Current function value: 0.574644923210144
    Current variable value: [ 0.82897025 -0.08847379  0.36032715  0.5334208 ] Current function value: 0.3903746008872986
    Current variable value: [ 0.7462303  -0.079105    0.3335184   0.47140023] Current function value: 0.2745198607444763
    Current variable value: [ 0.6791866  -0.07157812  0.31061673  0.42227766] Current function value: 0.19912949204444885
    Current variable value: [ 0.6240316  -0.06543481  0.29089484  0.38271144] Current function value: 0.14848771691322327


# Gradient descent using PyTorch's optmization


```python
x = Variable(torch.tensor([3.0,-1.0,0.0,1.0]),requires_grad=True)
  
optimizer = torch.optim.SGD([x], lr=1e-3, momentum=0.9) # create an optimizer that will do gradient descent optimization

for i in range(100):
  f = (x[0]+10.0*x[1])**2 + 5.0*(x[2]-x[3])**2 + (x[1]+2.0*x[2])**4 + 10.0*(x[0]-x[3])**4
  optimizer.zero_grad()
  f.backward()
  optimizer.step()
  if i%10==0:
    print("Current variable value:",x.detach().numpy(),"Current function value:", f.item())

```

    Current variable value: [ 2.694 -0.856  0.018  1.31 ] Current function value: 215.0
    Current variable value: [ 1.549712   -0.2964205   0.73650396  1.6820586 ] Current function value: 11.441824913024902
    Current variable value: [ 1.408362   -0.00583249  0.57971793  1.0532441 ] Current function value: 6.6437225341796875
    Current variable value: [ 0.6763986  -0.03913701  0.24533094  1.1681743 ] Current function value: 5.709615707397461
    Current variable value: [ 0.5996511  -0.13040473  0.5036107   0.53555894] Current function value: 1.093517541885376
    Current variable value: [ 0.59016085 -0.07009609  0.31826916  0.27739534] Current function value: 0.27947038412094116
    Current variable value: [ 0.5269876  -0.03624894  0.12591662  0.26216415] Current function value: 0.16900770366191864
    Current variable value: [ 0.45544916 -0.04479674  0.13050896  0.22580846] Current function value: 0.08700942993164062
    Current variable value: [ 0.399531   -0.04679754  0.18029977  0.18939461] Current function value: 0.034884385764598846
    Current variable value: [ 0.35905084 -0.03668797  0.18895635  0.19141704] Current function value: 0.022943012416362762


# Chain rule of derivative


```python
z = Variable(torch.tensor([1.0,-1.0]),requires_grad=True)

print("Variable z:",z.data)

def compute_x(z):
  x = torch.zeros(4)
  x[0] = z[0] - z[1]
  x[1] = z[0]**2
  x[2] = z[1]**2
  x[3] = z[0]**2+z[0]*z[1]
  return x

x = compute_x(z)
print("function x:",x.data)

def compute_f(x):
  f = (x[0]+10.0*x[1])**2 + 5.0*(x[2]-x[3])**2 + (x[1]+2.0*x[2])**4 + 10.0*(x[0]-x[3])**4
  return f

f = compute_f(x)
print("function f:",f.item())
print("")
# 
# Let's compute gradient of f with respect to x
g_x = torch.autograd.grad(f,x,retain_graph=True,create_graph=True)
print("Gradient of f with respect x:",g_x[0].data)
print("")
# Now compute Jacobian of x with respect to z and multiply with g_x to use chain rule
g_z = torch.autograd.grad(x,z,g_x,retain_graph=True) 

# But PyTorch can compute derivative of f with respect to z directly - this is the amazing capability!
g = torch.autograd.grad(f,z)

print("Gradient of f with respec to z by two-step calculation:",g_z[0].data)
print("Gradient of f with respec to z by one-step calculation:",g[0].data)
```

    Variable z: tensor([ 1., -1.])
    function x: tensor([2., 1., 1., 0.])
    function f: 390.0
    
    Gradient of f with respect x: tensor([ 344.,  348.,  226., -330.])
    
    Gradient of f with respec to z by two-step calculation: tensor([  710., -1126.])
    Gradient of f with respec to z by one-step calculation: tensor([  710., -1126.])


## Let's verify the result using math formula

$z_1 = 1, z_2 = -1$

$x_1 = z_1 - z_2 = 2, x_2 = z_1^2 = 1, x_3 = z_2^2 = 1, x_4 = z_1^2+z_1z_2 = 0$

## Let's compute partial derivative of $f$ with respect to $x_i$:

$\frac{\partial f}{\partial x_1} = 2(x_1+10x_2)+40(x_1-x_4)^3 = 2(2+10)+40(2-0)^3 = 24+320 = 344$

$\frac{\partial f}{\partial x_2} = 20(x_1+10x_2)+4(x_2+2x_3)^3 = 20(2+10)+4(1+2)^3 = 348$

$\frac{\partial f}{\partial x_3} = 10(x_3-x_4)+8(x_2+2x_3)^3 = 10(1)+8(3)^3 = 10+216 = 226$

$\frac{\partial f}{\partial x_4} = -10(x_3-x_4)-40(x_1-x_4)^3 = -10(1)-40(2)^3 = -10-320 = -330$

## Let's now compute partial derivative of $x_i$ with respect to $z_1$:

$\frac{\partial x_1}{\partial z_1} = 1$

$\frac{\partial x_2}{\partial z_1} = 2z_1 = 2(1) = 2$

$\frac{\partial x_3}{\partial z_1} = 0$

$\frac{\partial x_4}{\partial z_1} = 2z_1 + z_2 = 2(1)-1 = 1$

## Let's also compute partial derivative of $x_i$ with respect to $z_2$:

$\frac{\partial x_1}{\partial z_2} = -1$

$\frac{\partial x_2}{\partial z_2} = 0$

$\frac{\partial x_3}{\partial z_2} = 2z_2 = -2$

$\frac{\partial x_4}{\partial z_2} = z_1 = 1$

## Now we can apply the chain rule of derivative:

$\frac{\partial f}{\partial z_1} = \sum_i \frac{\partial f}{\partial x_i}\frac{\partial x_i}{\partial z_1} = 344(1) + 348(2) + 226(0) -330(1) = 710$

$\frac{\partial f}{\partial z_2} = \sum_i \frac{\partial f}{\partial x_i}\frac{\partial x_i}{\partial z_2} = 344(-1) + 348(0) + 226(-2) -330(1) = -1126$

## Using Jacobian notation, chain can be written as follows:

\begin{equation}
\nabla_z f = (J_z^x) (\nabla_x f) = 
    \begin{bmatrix}
        \frac{\partial x_1}{\partial z_1} & \frac{\partial x_2}{\partial z_1} & \frac{\partial x_3}{\partial z_1} & \frac{\partial x_4}{\partial z_1}\\
        \frac{\partial x_1}{\partial z_2} & \frac{\partial x_2}{\partial z_2} & \frac{\partial x_3}{\partial z_2} & \frac{\partial x_4}{\partial z_2}
    \end{bmatrix}
    \begin{bmatrix}
        \frac{\partial f}{\partial x_1}\\
        \frac{\partial f}{\partial x_2}\\
        \frac{\partial f}{\partial x_3}\\
        \frac{\partial f}{\partial x_4}
    \end{bmatrix} =
    \begin{bmatrix}
        1 & 2 & 0 & 1\\
        -1 & 0 & -2 & 1
    \end{bmatrix}    
    \begin{bmatrix}
        344\\
        348\\
        226\\
        -330
    \end{bmatrix} = 
    \begin{bmatrix}
        710\\
        -1126
    \end{bmatrix}
\end{equation}


# Optimization by gradient descent


```python
steplength = 1e-3 # for gradient descent
for i in range(1000):
  # function
  f = compute_f(compute_x(z))
  # Compute gradient of f with respect to z directly
  # PyTorch takes care of chain rule of derivatives
  g = torch.autograd.grad(f,z) 
  # adjust variable
  z = z - steplength*g[0]
  if i%100==0:
    print("Current variable value:",z.detach().numpy(),"Current function value:", f.item())
```

    Current variable value: [0.28999996 0.12600005] Current function value: 390.0
    Current variable value: [0.10041686 0.16472499] Current function value: 0.002078988589346409
    Current variable value: [0.0918402  0.16432461] Current function value: 0.0010676763486117125
    Current variable value: [0.08976527 0.1616178 ] Current function value: 0.0009482738096266985
    Current variable value: [0.08847897 0.15887022] Current function value: 0.0008560288697481155
    Current variable value: [0.08736441 0.15630378] Current function value: 0.0007776079582981765
    Current variable value: [0.08633856 0.15392138] Current function value: 0.000710221123881638
    Current variable value: [0.08538311 0.15170398] Current function value: 0.0006518377340398729
    Current variable value: [0.08448897 0.14963281] Current function value: 0.0006008768104948103
    Current variable value: [0.08364932 0.14769198] Current function value: 0.000556100916583091


# And of course optimization using PyTorch's gradient descent optimization


```python
z = Variable(torch.tensor([1.0,-1.0]),requires_grad=True)
  
optimizer = torch.optim.SGD([z], lr=1e-4, momentum=0.9) # create an optimizer that will do gradient descent optimization

for i in range(100):
  f = compute_f(compute_x(z))
  optimizer.zero_grad()
  f.backward()
  optimizer.step()
  if i%10==0:
    print("Current variable value:",z.detach().numpy(),"Current function value:", f.item())

```

    Current variable value: [ 0.929      -0.88740003] Current function value: 390.0
    Current variable value: [-0.19928937  0.4960087 ] Current function value: 0.9813485741615295
    Current variable value: [-0.4438892  0.7940718] Current function value: 24.753379821777344
    Current variable value: [-0.15691993  0.49109557] Current function value: 2.473496198654175
    Current variable value: [0.00094474 0.29169223] Current function value: 0.23995013535022736
    Current variable value: [0.0661319  0.20713682] Current function value: 0.024552151560783386
    Current variable value: [0.0927659  0.17456606] Current function value: 0.001717826584354043
    Current variable value: [0.10217324 0.16280487] Current function value: 0.0023682829923927784
    Current variable value: [0.10384294 0.1590711 ] Current function value: 0.003257464850321412
    Current variable value: [0.10240355 0.1583022 ] Current function value: 0.0029500178061425686


# Toy optimization example involving softmax and Huber loss function


```python
# A loss function measures distance between a predicted and a target tensor
# An implementation of Huber loss function is given below
# We will make use of this loss function in gradient descent optimization

# https://en.wikipedia.org/wiki/Huber_loss
    
def Huber_Loss(input,delta):
  m = (torch.abs(input)<=delta).detach().float()
  output = torch.sum(0.5*m*input**2 + delta*(1.0-m)*(torch.abs(input)-0.5*delta))
  return output

# Test Huber loss
a = torch.tensor([[0.3, 2.0, -3.1],[0.5, 9.2, 0.1]])
print(a.numpy())
ha = Huber_Loss(a,1.0)
print(ha.numpy())

b = torch.tensor([0.3, 2.0])
print(b.numpy())
hb = Huber_Loss(b,1.0)
print(hb.numpy())
```

    [[ 0.3  2.  -3.1]
     [ 0.5  9.2  0.1]]
    12.974999
    [0.3 2. ]
    1.545


# Here's the toy optimization problem

$min_{z_1,...,z_{10}} \{Huber\_Loss(p(z_1,...,z_{10})-y)\}$,

## where $y$ is a given 1-hot vector and each element of the vector $p$ is defined as:

$p_i = \frac{exp(z_i)}{exp(z_1)+...+exp(z_{10})}$

## The above function is known as softmax function. Loosely speaking, it turns any vector into  a probability vector. 


```python
y = torch.zeros(10)
y[2] = 1.0
print("Target 1-hot vector:",y.numpy())
z = Variable(torch.randn(y.shape),requires_grad=True) # Note that we are randomly initializing the z vector

optimizer = torch.optim.SGD([z], lr=1e-1, momentum=0.9) # create an optimizer that will do gradient descent optimization

# gradient descent
def gradient_descent(z,optimizer,softmax,loss,y,nIter,nPrint):
  for i in range(nIter):
    p = softmax(z)
    f = loss(p-y,1.0)
    optimizer.zero_grad()
    f.backward()
    optimizer.step()
    if i%nPrint==0:
      with np.printoptions(precision=3, suppress=True):
        print("Iteration:",i,"Variable:", z.detach().numpy(),"Loss: %0.6f" % f.item())
  return p,z
p,z = gradient_descent(z,optimizer,F.softmax,Huber_Loss,y,1000,100)
print("Optimum z:",z.data)
print("p at optimum:",p.data)
```

    Target 1-hot vector: [0. 0. 1. 0. 0. 0. 0. 0. 0. 0.]
    Iteration: 0 Variable: [ 1.727 -0.25   0.589  0.63   0.486  0.547  0.954  0.908  0.835  1.063] Loss: 0.488091
    Iteration: 100 Variable: [ 0.444 -0.381  5.313  0.244  0.163  0.199  0.384  0.368  0.34   0.416] Loss: 0.001682
    Iteration: 200 Variable: [ 0.41  -0.395  5.568  0.216  0.138  0.173  0.351  0.336  0.309  0.382] Loss: 0.000975
    Iteration: 300 Variable: [ 0.388 -0.404  5.73   0.199  0.122  0.156  0.331  0.316  0.29   0.361] Loss: 0.000688
    Iteration: 400 Variable: [ 0.372 -0.411  5.849  0.186  0.11   0.143  0.316  0.302  0.276  0.346] Loss: 0.000532
    Iteration: 500 Variable: [ 0.359 -0.416  5.944  0.176  0.1    0.134  0.304  0.29   0.264  0.333] Loss: 0.000433
    Iteration: 600 Variable: [ 0.349 -0.421  6.022  0.167  0.092  0.126  0.295  0.28   0.255  0.323] Loss: 0.000365
    Iteration: 700 Variable: [ 0.34  -0.425  6.089  0.16   0.086  0.119  0.286  0.272  0.247  0.315] Loss: 0.000315


    /usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:11: UserWarning: Implicit dimension choice for softmax has been deprecated. Change the call to include dim=X as an argument.
      # This is added back by InteractiveShellApp.init_path()


    Iteration: 800 Variable: [ 0.332 -0.428  6.148  0.154  0.08   0.113  0.279  0.265  0.24   0.307] Loss: 0.000278
    Iteration: 900 Variable: [ 0.326 -0.431  6.199  0.148  0.075  0.107  0.272  0.258  0.234  0.3  ] Loss: 0.000248
    Optimum z: tensor([ 0.3194, -0.4336,  6.2456,  0.1429,  0.0700,  0.1023,  0.2666,  0.2528,
             0.2282,  0.2946])
    p at optimum: tensor([0.0026, 0.0012, 0.9799, 0.0022, 0.0020, 0.0021, 0.0025, 0.0024, 0.0024,
            0.0026])


# Hessian computation (optional)


```python
x = Variable(torch.tensor([3.0,-1.0,0.0,1.0]),requires_grad=True)
f = (x[0]+10.0*x[1])**2 + 5.0*(x[2]-x[3])**2 + (x[1]-2.0*x[2])**4 + 10.0*(x[0]-x[3])**4

# PyTorch does not compute Hessian (second order derivatives) directly
# PyTorch can compute Jacobian vector product 
# We can use Jacobian vector product to compute Hessian

# Step 1: compute gradient
g = torch.autograd.grad(f,x,retain_graph=True,create_graph=True) # compute gradient with two important flags on

eye = torch.eye(4)

# Step 2: Use product of Jacobian of g and columns of identity matrix to compute hessian of f
H = torch.stack([torch.autograd.grad(g,x,eye[:,i],retain_graph=True)[0] for i in range(4)]) # hessian

print("Hessian:",H.data)
```

    Hessian: tensor([[ 482.,   20.,    0., -480.],
            [  20.,  212.,  -24.,    0.],
            [   0.,  -24.,   58.,  -10.],
            [-480.,    0.,  -10.,  490.]])


# Newton's method (optional)


```python
# Newton's optimization for an example - Powell Function (https://www.cs.ccu.edu.tw/~wtchu/courses/2014s_OPT/Lectures/Chapter%209%20Newton%27s%20Method.pdf)
# Minimize Powell function: f(x1,x2,x3,x4) = (x1+10x2)^2 + 5(x3-x4)^2 + (x2-2x3)^4 + 10(x1-x4)^4

x = Variable(torch.tensor([3.0,-1.0,0.0,1.0]),requires_grad=True)
f = (x[0]+10.0*x[1])**2 + 5.0*(x[2]-x[3])**2 + (x[1]-2.0*x[2])**4 + 10.0*(x[0]-x[3])**4

def LineSearch(x,d):
  minstep = 1.0
  minval=1e10
  for i in range(10):
    step = (i+1)/10.0
    xp = x.data + step*d.data
    fval = (xp[0]+10.0*xp[1])**2 + 5.0*(xp[2]-xp[3])**2 + (xp[1]-2.0*xp[2])**4 + 10.0*(xp[0]-xp[3])**4
    if fval < minval:
      minval = fval
      minstep = step
  return minstep

eye = torch.eye(4)

for itr in range(10):
  # Step 1: compute Newton direction d
  g = torch.autograd.grad(f,x,retain_graph=True,create_graph=True) # gradient
  H = torch.stack([torch.autograd.grad(g,x,eye[:,i],retain_graph=True)[0] for i in range(4)]) # hessian
  d = torch.solve(-g[0].unsqueeze(1), H)[0].t().squeeze() # solve Newton system
  
  # Step 2: update x with Newton direction d
  step_length = LineSearch(x,d)
  x.data += step_length*d.data # often step_length is set as 1.0
  print(x.data)
  f = (x[0]+10.0*x[1])**2 + 5.0*(x[2]-x[3])**2 + (x[1]-2.0*x[2])**4 + 10.0*(x[0]-x[3])**4
  print(f.data)
```

    tensor([ 1.5873, -0.1587,  0.2540,  0.2540])
    tensor(31.8025)
    tensor([ 1.0582, -0.1058,  0.1693,  0.1693])
    tensor(6.2820)
    tensor([ 0.7055, -0.0705,  0.1129,  0.1129])
    tensor(1.2409)
    tensor([ 0.4703, -0.0470,  0.0752,  0.0752])
    tensor(0.2451)
    tensor([ 0.3135, -0.0314,  0.0502,  0.0502])
    tensor(0.0484)
    tensor([ 0.2090, -0.0209,  0.0334,  0.0334])
    tensor(0.0096)
    tensor([ 0.1394, -0.0139,  0.0223,  0.0223])
    tensor(0.0019)
    tensor([ 0.0929, -0.0093,  0.0149,  0.0149])
    tensor(0.0004)
    tensor([ 0.0619, -0.0062,  0.0099,  0.0099])
    tensor(7.3712e-05)
    tensor([ 0.0413, -0.0041,  0.0066,  0.0066])
    tensor(1.4560e-05)


# Conjugate gradient method (optional)


```python
# Newton method using conjugate gradient (https://en.wikipedia.org/wiki/Conjugate_gradient_method)
# Hessian-vector product computation needs one autograd call per iteration inside CG iteration
# So this method never computes and stores the full hessian matrix
# CG solver might converge faster than other general linear equation solver
# Also, it seems to be more stable

x = Variable(torch.tensor([3.0,-1.0,0.0,1.0]),requires_grad=True)
f = (x[0]+10.0*x[1])**2 + 5.0*(x[2]-x[3])**2 + (x[1]-2.0*x[2])**4 + 10.0*(x[0]-x[3])**4

for itr in range(10):
  # Step 1: compute Newton direction d by CG method
  g = torch.autograd.grad(f,x,retain_graph=True,create_graph=True) # gradient
  r = -g[0].data.clone()
  p = r.clone()
  d = torch.tensor([0.,0.,0.,0.])
  rsold = torch.sum(r**2)
  for cg_itr in range(6): # cg_itr should be slightly larger length of variable - here variable length is 4
    q = torch.autograd.grad(g,x,p,retain_graph=True)[0] # hessian-vector (Jacobian-vector) product computation by autograd
    alpha = rsold/torch.sum(p*q)
    d += alpha*p
    r += -alpha*q
    rsnew = torch.sum(r**2)
    if rsnew<1e-10:
      break
    p = r + (rsnew/rsold)*p
    rsold = rsnew
    
  # Step 2: update x with Newton direction d
  step_length = LineSearch(x,d)
  x.data += step_length*d.data # often step_length is set as 1.0
  print(x.data)
  f = (x[0]+10.0*x[1])**2 + 5.0*(x[2]-x[3])**2 + (x[1]-2.0*x[2])**4 + 10.0*(x[0]-x[3])**4
  print(f.data)
```

    tensor([ 1.5873, -0.1587,  0.2540,  0.2540])
    tensor(31.8025)
    tensor([ 1.0582, -0.1058,  0.1693,  0.1693])
    tensor(6.2820)
    tensor([ 0.7055, -0.0705,  0.1129,  0.1129])
    tensor(1.2409)
    tensor([ 0.4703, -0.0470,  0.0752,  0.0752])
    tensor(0.2451)
    tensor([ 0.3135, -0.0314,  0.0502,  0.0502])
    tensor(0.0484)
    tensor([ 0.2090, -0.0209,  0.0334,  0.0334])
    tensor(0.0096)
    tensor([ 0.1394, -0.0139,  0.0223,  0.0223])
    tensor(0.0019)
    tensor([ 0.0929, -0.0093,  0.0149,  0.0149])
    tensor(0.0004)
    tensor([ 0.0619, -0.0062,  0.0099,  0.0099])
    tensor(7.3712e-05)
    tensor([ 0.0413, -0.0041,  0.0066,  0.0066])
    tensor(1.4561e-05)



```python

```

14aeedfbc4c199c7169672db4499ff35ae6da2a7
ipynb
Jupyter Notebook
Week 2/Partial_Derivatives_and_GD_and_chain_rule_new.ipynb
YongqinJian/Cmput-328
69f2b78f6c6b10d4445d5ba6bbff63bce0582b88
MIT
Week 2/Partial_Derivatives_and_GD_and_chain_rule_new.ipynb
YongqinJian/Cmput-328
69f2b78f6c6b10d4445d5ba6bbff63bce0582b88
MIT
Week 2/Partial_Derivatives_and_GD_and_chain_rule_new.ipynb
YongqinJian/Cmput-328
69f2b78f6c6b10d4445d5ba6bbff63bce0582b88
MIT
Qwen/Qwen-72B
1. YES
2. YES