In [1]:
import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

The data

In [2]:
x=np.array([-5, -4, -3, -2, -1,  0,  1,  2,  3,  4,  5])
y=np.array([-53.9, -28.5, -20.7,  -3.6,  -9.8,   5.0,   4.2,   5.1,  11.4, 27.4,  44.0])

Creating a matix with the powers of x

In [3]:
xmat=np.vstack((x**3,x**2,x,np.ones(len(x))))
xmat
Out[3]:
array([[-125.,  -64.,  -27.,   -8.,   -1.,    0.,    1.,    8.,   27.,
          64.,  125.],
       [  25.,   16.,    9.,    4.,    1.,    0.,    1.,    4.,    9.,
          16.,   25.],
       [  -5.,   -4.,   -3.,   -2.,   -1.,    0.,    1.,    2.,    3.,
           4.,    5.],
       [   1.,    1.,    1.,    1.,    1.,    1.,    1.,    1.,    1.,
           1.,    1.]])

The error function

In [4]:
def residuo(w):
    return (y-w.dot(xmat))
In [5]:
def err(w):
    return np.sum(residuo(w)**2)

Testing the function

In [6]:
err(np.array([0,0,0,0]))
Out[6]:
7140.3199999999997
In [7]:
np.sum(y**2)
Out[7]:
7140.3199999999997

ok

In [8]:
def derr(w):
    return -2.0*np.sum(residuo(w)*xmat,axis=1)

1. Gradient descent with epsilon=1.e-5

In [9]:
w=np.array([0,0,0,0])
print("erro: "+str(err(w)))
for i in range(50):
    w=w-1.0e-5*derr(w)
    if i%10==0:
        print("erro: "+str(err(w)))
print("erro: "+str(err(w)))
print("w: "+str(w))

erro: 7140.32
erro: 499.296253108
erro: 251.14181125
erro: 236.34401971
erro: 229.353674738
erro: 225.868405327
erro: 224.107804051
w: [ 0.40582178 -0.15346867  0.05946159 -0.00842238]
In [10]:
ypred= w[0]*x**3+w[1]*x**2+w[2]*x+w[3]
plt.plot(x, y, 'ro', label='Original data')
plt.plot(x, ypred, label='Fitted line')
plt.show()
In [11]:
# 2. Gradient descent with epsilon=1.e-4
In [12]:
w=np.array([0,0,0,0])
print("erro: "+str(err(w)))
for i in range(50):
    w=w-1.0e-4*derr(w)
    if i%10==0:
        print("erro: "+str(err(w)))
print("erro: "+str(err(w)))
print("w: "+str(w))

erro: 7140.32
erro: 357904.371903
erro: 5.36853782011e+22
erro: 8.0583576608e+39
erro: 1.20958686267e+57
erro: 1.81563097584e+74
erro: 5.22130488436e+89
w: [ -3.55918263e+42   7.75259950e+24  -1.69916825e+41  -2.16654622e+23]
In [13]:
ypred= w[0]*x**3+w[1]*x**2+w[2]*x+w[3]
plt.plot(x, y, 'ro', label='Original data')
plt.plot(x, ypred, label='Fitted line')
plt.show()

3. BFGS

In [14]:
sol1=minimize(x0=np.array([0,0,0,0]),fun=err,method="BFGS")
In [15]:
sol1
Out[15]:
      fun: 127.56506993008438
 hess_inv: array([[  8.09374160e-05,  -6.63256005e-10,  -1.44068549e-03,
          1.18721958e-08],
       [ -6.63256005e-10,   5.82750352e-04,   1.38931724e-08,
         -5.82750147e-03],
       [ -1.44068549e-03,   1.38931724e-08,   3.01896456e-02,
         -2.48662581e-07],
       [  1.18721958e-08,  -5.82750147e-03,  -2.48662581e-07,
          1.03729522e-01]])
      jac: array([  9.53674316e-07,   9.53674316e-07,   0.00000000e+00,
         0.00000000e+00])
  message: 'Optimization terminated successfully.'
     nfev: 66
      nit: 7
     njev: 11
   status: 0
  success: True
        x: array([ 0.29123927, -0.1799534 ,  2.45957738,  0.03589765])
In [16]:
ypred= sol1.x[0]*x**3+sol1.x[1]*x**2+sol1.x[2]*x+sol1.x[3]
plt.plot(x, y, 'ro', label='Original data')
plt.plot(x, ypred, label='Fitted line')
plt.show()

4. BFGS with my Jacobian

In [17]:
sol2=minimize(x0=np.array([0,0,0,0]),fun=err,method="BFGS",jac=derr)
sol2
Out[17]:
      fun: 127.56506993006992
 hess_inv: array([[  8.09375809e-05,  -4.39101880e-18,  -1.44068894e-03,
          5.25838054e-18],
       [ -4.39101880e-18,   5.82750583e-04,   9.22656049e-17,
         -5.82750583e-03],
       [ -1.44068894e-03,   9.23740251e-17,   3.01897177e-02,
         -1.26634814e-16],
       [  5.31259065e-18,  -5.82750583e-03,  -1.26634814e-16,
          1.03729604e-01]])
      jac: array([  2.84217094e-12,  -2.20268248e-13,   1.38555833e-13,
        -9.76996262e-15])
  message: 'Optimization terminated successfully.'
     nfev: 11
      nit: 7
     njev: 11
   status: 0
  success: True
        x: array([ 0.29123932, -0.17995338,  2.45957653,  0.03589744])

5. Nelder Mead

In [18]:
sol3=minimize(x0=np.array([0,0,0,0]),fun=err,method="Nelder-Mead")
sol3
Out[18]:
 final_simplex: (array([[ 0.29124112, -0.17995509,  2.45954951,  0.03587784],
       [ 0.29123943, -0.17995474,  2.45954791,  0.03590425],
       [ 0.29124162, -0.17994445,  2.45954219,  0.03579674],
       [ 0.29124209, -0.17994805,  2.4595116 ,  0.03588207],
       [ 0.29123556, -0.17994997,  2.45964811,  0.03591491]]), array([ 127.56506997,  127.56507001,  127.56507004,  127.56507004,
        127.56507006]))
           fun: 127.5650699702508
       message: 'Optimization terminated successfully.'
          nfev: 602
           nit: 360
        status: 0
       success: True
             x: array([ 0.29124112, -0.17995509,  2.45954951,  0.03587784])
In [ ]:
In [ ]: