title: gradient descent newton method using Hessian Matrix
body: <p>I am implementing gradient descent for regression using newtons method as explained in the 8.3 section of the Machine Learning A Probabilistic Perspective (Murphy) book. I am working with two dimensional data in this implementation. I am using following notations.<br>
      x =  input data points m*2<br>
      y =  labelled outputs(m) corresponding to input data<br>
      H = Hessian matrix  is defined as<br>
      </p>
  
  <p>gradient descent update<br>
      <img src="https://latex.codecogs.com/gif.latex?$$&space;%5Ctheta&space;=&space;%5Ctheta&space;-&space;H%5E%7B-1%7Dg(%5Ctheta)&space;$$" title="$$ \theta = \theta - H^{-1}g(\theta) $$" /></p>
  
  <p>where  loss function is defined as
      <img src="https://latex.codecogs.com/gif.latex?$&space;g(%5Ctheta)&space;=&space;2*&space;%5Csum_n&space;(%5Ctheta%5Et&space;x_%7Bn%7D&space;-&space;y_%7Bn%7D)x&space;$" title="$ g(\theta) = 2* \sum_n (\theta^t x_{n} - y_{n})x $" /></p>
  
  <p>In my case 
      <img src="https://latex.codecogs.com/gif.latex?$%5Ctheta$" title="$\theta$" /> is <img src="https://latex.codecogs.com/gif.latex?$2*1$" title="$2*1$" /> array and  H is <img src="https://latex.codecogs.com/gif.latex?$2*2$" title="$2*2$" /></p>
  
  <p>Here is my python implementation. However this is not working as cost is increasing in each iteration.</p>
  
  <pre><code>def loss(x,y,theta):
     m,n = np.shape(x)
     cost_list = []
     for i in xrange(0,n): 
       x_0 = x[:,i].reshape((m,1)) 
       predicted = np.dot(x_0, theta[i])    
       error = predicted - y
       cost = np.sum(error ** 2) /  m
       cost_list.append(cost)
  
     cost_list = np.array(cost_list).reshape((2,1))
     return cost_list
  
  
  def NewtonMethod(x,y,theta,maxIterations):
     m,n = np.shape(x)
     xTrans  = x.transpose()
     H       = 2 * np.dot(xTrans,x) / m
     Hinv    = np.linalg.inv(H)
     thetaPrev = np.zeros_like(theta)
     best_iter = maxIterations
     for i in range(0,maxIterations):
       cost = loss(x,y,theta)
       theta  = theta - np.dot(Hinv,cost))
       if(np.allclose(theta,thetaPrev,rtol=0.001,atol=0.001)):
          break;
       else:
         thetaPrev = theta
         best_iter = i
  
     return theta
  </code></pre>
  
  <p>Here are the sample values I used</p>
  
  <pre><code>import numpy as np
  
  x = np.array([[-1.7, -1.5],[-1.0 , -0.3],[ 1.7 ,  1.5],[-1.2, -0.7 ][  0.6,  0.1]])  
  y = np.array([ 0.3 ,  0.07, -0.2,  0.07,  0.03 ])
  theta = np.zeros(2)
  NewtonMethod(x,y,theta,100)
  </code></pre>
  
  <p>Need help / suggestions to fix this problem.<br>
  Thanks </p>