Implementing Logistic Regression with Scipy: Why does this Scipy optimization return all zeros?
Question: I am trying to implement a one versus many logistic regression as in Andrew
Ng's [machine learning class](https://www.coursera.org/learn/machine-
learning/home/week/5), He uses an octave function called `fmincg` in his
implementation. I have tried to use several functions in the
[`scipy.optimize.minimize`](http://docs.scipy.org/doc/scipy-0.14.0/reference/optimize.html),
but I keep getting all zeros in the classifier output, no matter what I put
in.

In the last many hours, I've checkout out a ton of resources, but the most
helpful have been [this stack overflow
post](http://stackoverflow.com/questions/18801002/fminunc-alternate-in-numpy),
and [this blog post](http://aimotion.blogspot.com/2011/11/machine-learning-
with-python-logistic.html).

Is there any obvious or not so obvious place where my implementation has gone
astray?

    
    
    import scipy.optimize as op
    
    def sigmoid(z):
        """takes matrix and returns result of passing through sigmoid function"""
        return 1.0 / (1.0 + np.exp(-z))
    
    def lrCostFunction(theta, X, y, lam=0):
        """
        evaluates logistic regression cost function:
    
        theta: coefficients. (n x 1 array)
        X: data matrix (m x n array)
        y: ground truth matrix (m x 1 array)
        lam: regularization constant
        """
        m = len(y)
        theta = theta.reshape((-1,1))
        theta = np.nan_to_num(theta)
    
        hypothesis = sigmoid(np.dot(X, theta))
        term1 = np.dot(-y.T,np.log(hypothesis))
        term2 = np.dot(1-y.T,np.log(1-hypothesis))
        J = (1/m) * term1 - term2
        J = J + (lam/(2*m))*np.sum(theta[1:]**2)
        return J
    
    def Gradient(theta, X, y, lam=0):
        m = len(y)
        theta = theta.reshape((-1,1))
        hypothesis = sigmoid(np.dot(X, theta))
        residuals = hypothesis - y 
        reg = theta  
        reg[0,:] = 0
        reg[1:,:] = (lam/m)*theta[1:]
        grad = (1.0/m)*np.dot(X.T,residuals)
        grad = grad + reg
        return grad.flatten()
    
    def trainOneVersusAll(X, y, labels, lam=0):
        """
        trains one vs all logistic regression.
    
        inputs:
            - X and Y should ndarrays with a row for each item in the training set
            - labels is a list of labels to generate probabilities for.
            - lam is a regularization constant
    
        outputs:
            - "all_theta", shape = (len(labels), n + 1) 
        """
        y = y.reshape((len(y), 1))
        m, n = np.shape(X)
        X = np.hstack((np.ones((m, 1)), X))
        all_theta = np.zeros((len(labels), n + 1))
        for i,c in enumerate(labels):
            initial_theta = np.zeros(n+1)
            result, _, _ = op.fmin_tnc(func=lrCostFunction, 
                                fprime=Gradient,
                                x0=initial_theta, 
                                args=(X, y==c, lam))
            print result
            all_theta[i,:] = result
        return all_theta
    
    
    def predictOneVsAll(all_theta, X):
        pass
    
    a = np.array([[ 5.,  5.,  6.],[ 6.,  0.,  8.],[ 1.,  1.,  1.], [ 6.,  1.,  9.]])
    k = np.array([1,1,0,0])
    # a = np.array([[1,0,1],[0,1,0]])
    # k = np.array([0,1])
    solution = np.linalg.lstsq(a,k)
    print 'x', solution[0]
    print 'resid', solution[1]
    thetas = trainOneVersusAll(a, k, np.unique(k))
    


Answer: The problem lies in your `Gradient` function. In `numpy` assignment is **not
copying objects** , so your line

    
    
    reg = theta
    

makes `reg` a reference to `theta`, so each time you compute gradient you
actually modify your current solution. It should be

    
    
    reg = theta.copy()
    

I would also suggest starting from random weights

    
    
    initial_theta = np.random.randn(n+1)
    

Now solution is no longer zeros (although I did not check each formula, so
there still might be mathematical error). It is also worth noting, that for
linearly separable problems, logistic regression without regularization is
ill-posed (its objective is unbounded), so I suggest testing with `lam>0`.

