Logistic regression and singular Hessian

Question

I've been following along with Andrew Ng's excellent course on Machine Learning (CS 229), and have been working on Problem Set 1.

I'm trying to fit a logistic regression model using Newton's Method. I have tried several approaches:

1) Calculating the components of $\mathbf{\theta} := \mathbf{\theta} - \mathbf{H}^{-1}\mathbf{\nabla}_{\mathbf{\theta}}\mathbf{\ell}$ element-by-element then solving;

2) Updating $\mathbf{\theta}$ using $\mathbf{(X^{T}WX)^{-1}X^{T}W\mathbf{z}}$ where $\mathbf{z}:=\mathbf{X}\theta+\mathbf{W^{-1}(y-p)}$. Described on slide 21 here.

3) Using the scikit's built-in package LogisticRegression to solve the system.

Number 1 gives me a singular Hessian. Number 2 gives a singluar weightings matrix $\mathbf{W}$ after 2 or 3 iterations. Number 3 converges to $\theta=[-1.0035,0.6191,1.0792]^{T}$. My data matrix $\mathbf{X}$ has 3 columns $[x_0,x_1,x_2]$, where $x_0$ is a column of ones. What am I doing wrong?

From a plot of the data (below), I don't believe the features are linearly dependent.

Code for approach number 2 below. The dataframe 'data' has four columns $[x_0,x_1,x_2,y]$.

            import numpy as np
            import numpy.linalg as linalg
            import pandas as pd
            import csv as csv

            # Read in data to a pandas dataFrame. This file has 3 columns, 'x1', 'x2' and 'y'
            data = pd.read_csv(r'.\logistic_xy.csv')

            # Convert to a numpy array
            data = np.array(data)

            # Define data matricies

            X = np.array(data[0::,0:3])
            y = np.array(data[0::,3])
            y = y.reshape((99,1))
            W = np.zeros([99,99],dtype=float)
            p = np.zeros([99,1],dtype=float)
            z = np.zeros([99,1], dtype=float)

            # Define hypothesis function
            def hTheta(x):
                g = np.exp(x)/(1+np.exp(x))
                return g

            # Define convergence variable and initialize theta

            theta = np.zeros([3,1],dtype=float)
            thetaNew = theta
            conv = 1

            while conv > 0.01:
                for i in xrange(0,int(data[0::,0].sum())):
                    theta = thetaNew
                    x = np.dot(np.transpose(theta),data[i,0:3])[0]
                    h = hTheta(x)
                    W[i,i] = h*(1-h)
                    p[i] = h

                z = np.add(np.dot(X,theta),np.dot(np.linalg.inv(W),np.subtract(y,p)))
                A = np.linalg.inv(np.dot(X.transpose(),np.dot(W,X)))
                B = np.dot(X.transpose(),np.dot(W,z))

                thetaNew = np.dot(A,B)
                conv = np.linalg.norm(np.subtract(thetaNew,theta))